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NUMERICAL SIMULATION OF A COMPRESSIBLE, 

HOMOGENEOUS, TURBULENT SHEAR FLOW 

Abstract 

A direct, low Reynolds number, numerical simulation has been per- 
formed on a homogeneous turbulent shear flow. The full compressible 
Navier-Stokes equations were used in a simulation on the ILLIAC IV com- 

3 

puter with a 64 mesh . 

The flow fields generated by the code were used as an experimental 
data base, to examine the behavior of the Reynolds stresses in this 
simple, compressible flow. Emphasis was placed on determining the 
variation of the structure of the stresses and their dynamic equations 
as the character of the flow changed . 

The objectives of this work center round the modeling of these 
stresses in a compressible turbulent flow. It has been found that the 
structure of the stress tensor is more heavily dependent on the shear 
number and less on the fluctuating Mach number than was originally 
thought . 

The pressure-strain correlation tensor in the dynamic stress equa- 
tions can be directly calculated in this simulation. It is found that 
these correlations can be decomposed into several parts, as contrasted 
with the traditional incompressible decomposition into two parts . The 
performance of existing models for the conventional terms is examined, 
and a new model is proposed for the "mean-fluctuating" part . 

The additional terms in the pressure-strain tensor relate to the 
compressibility of the fluid. They are found to be of the same order of 
magnitude as the conventional terms. The modeling, therefore, becomes 
quite Important in an averaged simulation of the Navier-Stokes equa- 
tions. The behavior of these terms is examined and suggestions made for 
their modeling. A new class of models based on a structural similarity 
concept is examined. This type of model is used for the entire pressure 
strain tensor and looks promising. 

In these simulations the computer is used as a numerical wind tun- 
nel, a capability only recently available with the advent of large. 


iv 



fast, "vector" machines like the ILLIAC IV. The ability to measure 
quantities that were previously Inaccessible from these simulations 
should prove to be a great boon to turbulence research. Further direct 
and large eddy simulations In "building block" flows , such as the homo- 
geneous shear flow, are therefore highly recommended. They can lead to 
a fuller understanding of the physics of turbulence, while at the same 
time providing Information to the conventional turbulence modeler who Is 
Interested In technologically useful flows . 
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Chapter I 


MOTIVATION AND LITERATURE 


1 .1 Introduction 

Modeling of turbulent fluid flows is a subject of great scientific 
and technological interest . Turbulent flows by their nature contain a 
large range of sizes of motion (length scales) and an equivalent range 
of time scales. This range of scales causes problems for the flow simu- 
lator because the smallest and largest scales of motion in his flow 
cannot be represented at the same time on a relatively coarse mesh. The 
required mesh systems are beyond the capability of modern computers. 

Ife believe that the Navier-Stokes equations describe the flow of a 
Newtonian fluid. Reynolds (1883) applied an averaging operator to the 
Navier-Stokes equations in the hope that the resulting equations would 
be easier to solve. Because these equations are nonlinear, averaging 
introduces unknown correlations that prevent the system of equations 
from being complete (closed), unless assumptions are introduced about 
hov? these correlations behave. The unknown correlations in the momentum 
equations are the Reynolds stresses and assumptions about their behavior 
are turbulence models . 

Turbulence modeling has received much attention over the years . 
Most of the fundamental work in this area has been based on the imcom- 
presslble Navier-Stokes equations. However, most flows of technological 
interest are compressible. For many years, models developed from the 
Incompressible equations have been applied to compressible flows, in 
many cases with great success. At higher Mach numbers, existing turbu- 
lence models become increasingly Inadequate. As it is very difficult to 
make experimental measurements in these flows we do not know the reason 
for this failure . 

It isi the purpose of this work to study what happens to the Rey- 
nolds stresses at high Mach number and investigate how they can be 
modeled. Wfe use direct simulation of the full, unaveraged, Navier- 
Stokes equations to study this problem. Even with the power of a modern 
vector computer we are limited to simple geometries and low Reynolds 
numbers. By simulating the full Navier-Stokes equations we have no 
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closure problem. Ife can introduce an averaging process on the simulated 
results and directly calculate the various turbulence terms. Simula- 
tions such as these give us insight and guidance for constructing turbu- 
lence models that will be of use in more complicated flows . 

In this work we perform simulations of compressible homogeneous 
turbulent shear flows on the ILLIAC IV computer and use the resulting 
flow fields as data bases . 

Our objectives are : 

1. To develop techniques for direct simulation of compressible 
homogeneous turbulent flows . 

2 . To apply these techniques to a shear flow simulation on the 
ILLIAC IV. 

3 . To develop a data base from these simulations . 

4. To study the Reynolds stresses and terms in the Reynolds 
stress equations searching for the effects of compressibil- 
ity . 

5. To test turbulence models by comparing them with exact 
results computed from the data base. 

In the remainder of this chapter we discuss the origion of the 
Reynolds stress, the various types of averaging and the complications 
that arise in compressible flows . Ivfe describe previous attempts to 
understand compressible turbulence and further discuss the reasons that 
led us to a direct simulation. 

In chapter II we present the mathematical foundations of these 
simulations and select numerical method . 

In chapter III we show how the equations and the numerical methods 
are implemented in a computer code on a vector computer. kfe describe 
the testing of the computer codes and the time development of a typical 
simulation. 

In chapter IV we present results. Vfe first describe character- 
ization of these simulations. Then we present the data base. Vte dis- 
cuss the structure of the Reynolds stresses and then some of the terms 
in their dynamic equations . Wfe evaluate several turbulence models and 
propose some improvements. 

Chapter V contains the conclusions . 
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1 .2 Averaging and the Origin of the Reynolds Stresses 


Solutions of the Navler-Stokes equations represent the flow of a 
Newtonian fluid. They are the governing equations for laminar and tur- 
bulent flows at all speeds. As alluded to earlier, structures in 
turbulent flows , at technologically interesting Reynolds numbers , con- 
tain too large a range of sizes for representation with current compu- 
ters, 

Hirt (1969) has shown that the number of mesh points in a simple, 
complete (no turbulence model) , three-dimensional simulation scales on a 
Reynolds number formed from the turbulence quantities, 

N^ « (1.2.1) 

where q is a turbulent vleocity scale, L is a turbulent length 

scale, and v is the kinematic viscosity. N is the requred number of 

O 

mesh points in one direction. For large N is beyond current 

storage capacity. This necessitates that a smoothing or averaging be 
applied to reduce this range in length scales so that we may stay within 
the resolution capacity of modern computers. This averaging is normally 
applied directly to the Navier-Stokes equations. Consequently, an 
averaged form of these equations, including a turbulence model, is 
usually solved . 

We write the full Navler-Stokes equations using tensor notation: 
Conservation of Mass 


P,t+(PUi),i = 0 (1.2.2) 

Conservation of Momentum 




Conservation of Energy 


(1.2.3) 


(pe) + (u. (pe+P)) . 
3 -*• > J 


^'"i'^ij\j" ^ j,j 


(1.2.4) 
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where 


Pe 


VFT) 


+ p 




= pc + 


PUiUi 


(1.2.5) 


is the sum of the Internal and kinetic energies per unit volume . Wfe 
have used K for the internal energy per unit mass because the standard 
symbol u is used for velocity. 


'ij - -T 


( 1 . 2 . 6 ) 


is the viscous stress tensor, and 




= - kT 


JJ 


(1 .2.7) 


is the heat flux vector. The pressure P is related to C and p 
through an equation of state, P = f(5,p)> in our case the perfect gas 
law with constant specific heats, kb divide the variables into mean and 
fluctuating components , 


p = p + p' ; u^ = Uj, + uj^ ; e = e + e' ; P = P + P' 

( 1 . 2 . 8 ) 


and apply Reynolds averaging to equations (1 .2 .2)-(l .2 .7) . The averag- 
ing operator is as yet unspecified, but has the property that ulj^ = 0. 
Writing the momentum equation only and indicating the average with a 
bar, 


(P + (P + 





|^p(u|up + u^(p'u^) + Ujp^uJ 
Reynolds Stress 


(1.2.9) 


we see the appearance of additional terms , the Reynolds stresses . These 
terms Increase the number of variables in the problem and require addi- 
tional equations if we are to have a complete set . Ihe term with the 
time-derivative now appears as two terms; the second term must also be 
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modeled. Seeking additional equations, we can form equations for the 
Reynolds stresses, but unfortunately we find the appearance of yet 
higher-order correlations in these equations. This process of forming 
equations for these higher-order correlations is hopeless, because we 
can never complete the set . Via must Introduce a model for the Reynolds 
stresses that relates them to the mean-flow variables . 


A standard simplification is to assume that the flow has constant 
density and is therefore incompressible (u. = 0) . The Reynolds 
stresses in Eq . (1.2.9) then reduce to u|uj ,and the second time- 
derivative term disappears. Favre (1965) introduced a variation of the 
averaging procedure that simplified the appearance of the Reynolds 
stresses in the compressible flow. He multiplied the velocity and 
energy by the density and averaged this product to form "mass-weighted" 
variables . 


Pu. 


u 1 “ '^i ~ ^1 


( 1 . 2 . 10 ) 


Tills averaging has the property pu^ = 0. Applying this definition to 
Eqs. (1.2.2) through (1.2.7), we find that the appearance of the Rey- 
nolds stresses has been simplified and looks quite similar to the incom- 
pressible stress. 


p ulu'. + u.p'ul + u.p'u'. - p'uju'. 
ij ij ij 

Reynolds-averaged 


pij^uT (1 .2 .11) 

Favre- averaged 


Thie simplification is in appearance alone and has nothing to do with the 
physics . It is the form that most modern compressible flow simulations 
use. Ife shall show that Reynolds- and Favre-averaging are identical in 
the homogeneous flows that we simulate, and that our conclusions apply 
to either type. 

Favre-averaging of the Navler-Stokes equations is thoroughly 
treated in Rubesin and Rose (1973), to which the reader is referred for 
a complete discussion. 
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1 .3 A Short History 


Early studies of compressible turbulence applied decompositions to 
the flow variables and studied their interaction. Kovasznay (1953) and 
Chu and Kovasznay (1958) proposed the idea of "modes." Ihey derived 
equations for the vorticity mode, corresponding to convective and vorti- 
cal motions, and for the pressure (acoustic) and entropy (temperature) 
modes. They analytically studied interactions among these modes, but 
were limited to low turbulence levels by their analytic techniques . 
They concluded that these interactions were second-order and therefore 
quite small at low- turbulence levels. Moyal (1952) concluded the same 
in his analysis. He divided the kinetic energy spectra into "eddy 
turbulence" and "random noise" parts. He did this by Fourier- 
transforming the velocity vector and decomposing it into vectors that 
are, respectively, perpendicular and parallel to the wave-number 
vector. His analysis, like that of Qiu and Kovasznay, was limited to 
low turbulence levels and therefore not applicable in technologically 
useful flows . 

We recall this work because we are now able to perform these decom- 
positions numerically, without the restriction to low turbulence levels. 
Wfe shall use Moyal' s decomposition in the presentation of some results 
In Chapter IV. 

Because measurements in high-speed flows are so difficult, little 
is known about the structure of the Reynolds stresses and their equa- 
tions. Historically, simulators have used models derived from the in- 
compressible Navier-Stokes equations, torkovin (1962) used the limited 
data available at the time to show that the Re3molds stresses in super- 
sonic boundary layers were structurally similar to incompressible flow. 
Laufer (1969) used Favre averages to come to the same conclusion. Over 
the years, application of incompressible models has met with great suc- 
cess in boundary layers. However, limits on this applicability began to 
be recognized. Bradshaw (1977) quantified Morkovin's hypothesis. He 
agreed with Morkovln that incompressible models should not be applied in 
boundary layers with external Mach numbers greater than five, nor in 
boundary layers with large pressure gradients (shock- boundary layer in- 
teraction) . Bradshaw also concluded that these models are inappropriate 
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in. free-stream layers with Mach numbers greater than 1.5; however, no 
physical reason Is given. A general conclusion of the Free Shear Layers 
Conference (Morkovin et al. (1972)) was that existing turbulence models 
were inadequate for these flows . None was able to predict the well- 
known decrease in mixing-layer spreading rate with increased Mach niomber 
without introducing empirical data. 

These failures have variously been attributed to variations in den- 
sity or some unspecified "compressibility effect ." Brown and Roshko 
(1974) performed a low-speed, variable density, mixing layer experiment. 
Their flow was virtually incompressible, but they controlled the density 
ratio of the two streams by using gases of different molecular weights . 
They found no evidence of spreading rate variation with density ratios 
appropriate for single-component , high-speed mixing layers and concluded 
that spreading-rate variation must be a "compressibility effect ." Oh 
(1974) simulated a two-dimensional, high-speed mixing layer. He pro- 
posed a kinetic-energy equation model that provided for a nonzero 
pressure-dilatation interaction. He was able to correctly predict the 
trend of the spreading, leading us to suspect that there are unrecog- 
nized physical phenomena that must be modeled . The questions that re- 
main unanswered for lack of experimental data are, " Wiat changes occur 
in the Reynolds stresses and their equations in a compressible flow, and 
how should they be modeled?" In the next section, we discuss how we 
approached these questions . 

1 .4 The i^proach of This kbrk 

Much progress in turbulence research has come from the study of 
homogeneous flows. These are flows that extend to infinity in all di- 
rections and are statistically similar everywhere in space. Presumably 
a homogeneous flow is an approximation to a piece of an inhomogeneous 
flow. It allows us to separate and distinguish competing processes in 
the development of the turbulence and is also amenable to analytic 
treatment (Batchelor, 1953). 

Several homogeneous turbulent-flow experiments have been performed 
(box-turbulence: Comte-Bellot and Corrsin, 1971, and Bennett, 1976; 
shear flow: Rose, 1966, Champagne et al., 1970, and Ifarris et al . , 
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1977); plane strain: llicker and Reynolds, 1975) to obtain experimental 
measurements of turbulence quantities. All of these experiments were 
performed in steady-flow wind tunnels and substituted spatial (down- 
stream) development for time development. Consequently, they are ap- 
proximations to homogeneous flow and are slightly Inhomogeneous in the 
downstream direction. Analysis by Harris (1977) showed that this shear 
flow has about the maximum mean-velocity gradient possible while still 
maintaining an approximation to homogeneity. Ihe downstream distance in 
which a significant change In turbulence quantities occurs must be 
larger than the largest turbulent length scales in the flow. This can- 
not be satisfied at high speeds with large velocity gradients. There- 
fore, there are do high Mach number homogeneous flow experiments. 
However, we would still like to study a compressible homogeneous flow. 
This is where the work must start if we are to understand turbulence in 
a compressible flow. 

In the last few years, the advent of very large computers has pro- 
vided the capability of doing three-dimensional simulations of flows 
with simple geometries (Deardorff, 1970; Orszag, 1971; Qark, 1977; 
Mansour, 1978; Moin, 1978; and Pulliam, 1979). The majority of turbu- 
lent flow simulations have solved averaged or filtered equations and 
incorporated a turbulence model. 

Because of the absence of experimental measurements of compressible 
turbulence, we use a large vector computer as a numerical wind tunnel 
and perform numerical experiments . This approach was pioneered by dark 
(1977). He simulated box turbulence with the unaveraged, incompressible 
Navier-Stokes equations and therefore did not need a model to close his 
equations. In contrast to experiment, his simulations (and all homogen- 
eous flow simulations) developed in time and was therefore a closer ap- 
proximation to the ideal. He compared the performance of turbulence 
models against the "exact" terms as calculated from his simulated flow 
fields. 

In this work, we simulate a compressible, homogeneous, turbulent 
shear flow by solving a transformed version of the full Navier-Stokes 
equations ((1.2.2) through (1.2.7)). A shear flow is the simplest 
turbulent flow with a continuous source of turbulent kinetic-energy 
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production. At low Mach number we may also compare with the shear flow 
experiments listed above. 

The mean flow Is pictured In Fig. 1.1. Geometrically, It Is like a 
deck of playing cards, with each card sliding on top of the one below 
It . 


1 .5 Numerical Requirements 

NS cannot simulate a homogeneous flow with an Infinite domain. 
Like an experimentalist, we choose a portion of the flow and Impose 
boundary conditions . NS must choose a large enough portion of the flow 
that the turbulent length scales are much smaller than our domain, or 
the boundary conditions will Interfere. Ideally, we would like to Im- 
pose periodic boundary Conditions, l.e., one side of our computational 
mesh Is Identical to the opposite side. These boundary conditions are 
essentially transparent to the simulation. If the turbulent length 
scales are small enough. However, Inspection of Fig. 1.1 shows that the 
shear flow Is not periodic. The mean velocity varies across a computa- 
tional mesh. In Chapter II we shall Introduce an analytic coordinate 
transformation on the Navler-Stokes equations that allows the equations 
to have periodic solutions . 

If we are to use the simulated flow fields as a data base, we must 
be sure of their accuracy. Nfe require a highly accurate numerical 
method In both space and time. This requirement Is also addressed In 
Chapter II, where we discuss the equations and the numerical method. 
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Chapter II 


EQUATIONS AND NUMERICAL METHOD 

In this chapter we discuss the coordinate transformation that we 
apply to the Navier-Stokes equations . Ife describe an alternate way of 
writing the momentum equations that ensures the numerical conservation 
of quantitie that are analytically conserved, and we justify the choice 
of a numerical method . 

2 .1 The linear Coordinate Transformation 

Wfe apply a coordinate transformation to the Navier-Stokes equations 
that allows them to have periodic solutions . fe then impose Periodic 
Boundary (bnditions (PBCs) on the transformed equations and use Fourier 
methods for the spatial derivatives . 

Homogeneous flows extend to infinity in space . Obviously we cannot 
simulate the entire flow, nor is it necessary. Vfe choose a portion of 
the flow field and impose boundary conditions on the edge of this do- 
main. Boundary conditions are a source of numerical uncertainty in any 
simulation. In simulations of a homogeneous flow, we may reduce this 
uncertainty by applying PBCs after a suitable coordinate transformation. 
PBCs enforce all variables to be periodic on the domain. They are es- 
sentially transparent to the simulation if the domain is significantly 
larger than the largest turbulent length scales that naturally develop 
inside it. PBCs are desirable because they do not introduce unknown 
effects into the flow. They are the least Intrusive boundary conditions 
available for this geometry, and we would like to use them. 

^fost homogeneous flows are not periodic (except isotropic flows) , 
because the mean velocity is not periodic. This can be seen by dividing 
the variables into mean and fluctuating components . kfe do this in gen- 
eral for all homogeneous flows with linear mean-velocity gradients, and 
then show how the derivation is made specific for homogeneous shear 
flows. Vfe write: 
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p' + Po + EjXj ; 


u. 


u! + A 


Ij j 


P = 


P' + P + F.x, 
o J j 
( 2 . 1 . 1 ) 


where the mean parts of the density and pressure have both constant and 
linear gradient parts. Tensor and vectors Ej and Fj represent 
the linear gradients of each variable. They are constants for our case 
but may be generalized to vary with time. For example, in the shear 
flow. 


"Ij 


auf 


"0 S 0 
0 0 0 
0 0 0 


( 2 . 1 . 2 ) 


where S is the "shear rate," 3U/3y. fe shall show that Ej and 
Fj must be zero when we introduce the coordinate transformation. 

If we introduce the definitions (2.1.1) into the conservation of 
mass equation (2.1.2), we immediately produce an equation with the coor- 
dinates as coefficients. 




(2.1.3) 


+u' (p +Ex)+A (p'+p +Ex)+Eu' = 0 

i,i o j il^ o j j' j j 


Equation 2.1.3 will not allow periodic solutions, because the coordi- 
nate, Xj, appears explicitly in the coefficients. Hence we are led to 
apply a coordinate transformation in order to eliminate these coeffici- 
ents and to derive equations for the fluctuating (primed) quantities 
only. 

Ife define this transformation for the general case as 


= %jXj, t’ = t 


(2.1.4) 


where relates the transformed corodlnates to the Cartesian coor- 
dinate. This idea originated with Rogallo (1979), who first applied it 
to the solution of an incompressible homogeneous shear flow. He substi- 
tuted the decomposed field (2.1.1) and the transformation (2.1.4) Into 
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the incompressible momentum equations and derived the conditions that 
allow the equations to have periodic solutions . Ivfe follow t^e same 
reasoning for the full Navier-Stokes equations . The same procedure is 
followed for all the equations . fe illustrate this for the mass equa- 
tion only because the algebra is very messy. 

Inserting (2.1.1) and (2.1.4) into the Navier-Stokes equations 
(1 .1 .2)-(l .1 .7) and writing the mass equation, the derivatives become 


3 „ 3 

3x^ ^ ^ 


3 3 • ,3 

Jt S'P' ®jk^ki^i IxT 

j 


where is the inverse of . 


(2.1.5) 




Sk^ ’ \ jk\^mk'^ j 


(Bjk + P' ■ \ 


®kj (^mi ^ '"j 


+ A, C ,C , + A. C -C , ) 
jn nx mk jn m£ nk/ 


jn 


X' 


( 2 . 1 . 6 ) 


Ivfe have isolated all the coordinate coefficients in the last two 
terms. Ideally, we would choose some Bjj^ and the last two terms in 
(3.1.6) would be zero. However, there is no Bjj^ that will zero the 
coefficient in the last term. Therefore, must be identically zero, 

i.e., there can be no mean gradient of density. 

The second-to-last term has a coefficient that we may set equal to 
zero, resulting in a set of coupled ordinary differential equations: 

Solution of these equations, subject to the initial conditions. 


®lj “ ^Ij (Cartesian mesh) 


( 2 . 1 . 8 ) 


defines the transform for a specific mean velocity gradient tensor. 
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As stated, the same procedure Is followed to transform the momentum 
equations. Equation (2.1.7) also results independently from the momen- 
tum equations . 

Mien applying the procedure to the total energy equation (1.2.4), 
we were unable to zero all of the coordinate coefficients. The reason 
can be seen by inspecting the definition of the total energy, Eq . 
(1.2.5), the total energy is the sum of the kinetic and internal ener- 
gies and is quadratic in (2.1.1) and the velocity. This leaves some 
terms in (1.2.4) quadratic and some cubic when Eqs. (2.1.5) are inser- 
ted. A linear coordinate transformation cannot eliminate the coordinate 
coefficients from all these terms at the same time . 

Vfe solve this problem by subtracting the kinetic energy equation 
from (1.2.4), leaving only the thermal energy equation. fe may do this 
and retain the complete set of equations, because the kinetic energy 
equation is not independent of the mass and momentum equations . 


Taking the time derivative of (1.2.5), we find 


3 3 P . 9 “i'^i 

3t “ 3t (Y-1) 3t ^ 2 


(2.1.9) 


which shows that the total energy equation is the svun of the internal 
energy P/(y-l) and the kinetic energy per unit volume. Ife subtract 
the kinetic energy equation and multiply by (y- 1) to form the equation 
for the pressure before we apply the transformation. 


3 

Tt 


P +YP„^_l 


+ u. P j 

1 .1 


(,-l) - (v-1) (2.1.10) 


The total energy is analytically conserved . By solving the thermal 
energy equation (2.1.10), we give up guaranteed numerical conservation 
of total energy, but we shall show later how this is regained. 

Equation (2.1.10) can now be transformed in the same manner as the 
mass and momentxim equations. Without presenting the complicated alge- 
bra, we find that Fj must be identically zero, for the same reason 
as found with the density. There can be no mean gradient in the pres- 
sure. This is in contrast to the incompressible case, where the mean 
pressure gradient is arbitrary. This is due to the fact that the 
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pressure is a true flow variable and not the result of a kinematic con- 
straint, as it is in incompressible flow. 


This derivation is general for an arbitrary mean velocity gradient, 
Aji j . Taking the corresponding to the shear flow, Eq . (2.1.2), and 


solving for the transform tensor, 
initial condition (2.1.8), we find 

~1 -St 




1 

0 


or 




X 


y' = y 


in (2.1.7) subject to the 


St y 




( 2 . 1 . 11 ) 


( 2 . 1 . 12 ) 


= z 


This is the coordinate system for the homogeneous shear flow that allows 
the equations for the fluctuating quantities to have periodic solutions . 
It is effectively "glued" on top of the mean velocity of Fig. 1.1, and 
shears along with it. 

Before presenting the final form of the equations that are used in 
the simulation, we show an alternative form of the momentum equations 
(1.2.3) that regains the total energy conservation lost by solving the 
thermal energy equation. This form will also guarantee that kinetic 
energy is not being artificially produced . 


2 .2 The Ojnservatlon Properties 

From previous experience (Msmsour, 1978), it is known that artl-^ 
ficlal generation of conserved flow-field quantities by the numerical 
method can destroy the validity of a simulation. For example, in in- 
compressible flows, in the absence of viscosity and turbulent kinetic 
energy production, and, with PBCs , it can be proven that the total mass, 
momentum, and kinetic energy remain constant. Numerical simulations do 
not always ensure this. Rwak (1975) and Shaanan (1975) used a modified 
(but exact) form of the convective terms in the momentum equations, in 
order to ensure kinetic energy conservation under these conditions. 
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In a compressible flow under the same conditions (periodic boundary 
conditions, no turbulent kinetic energy production, zero viscosity), 
kinetic energy conservation does not hold . There is a physical mecha- 
nism for exchange of energy between the kinetic energy and the thermal 
energy through a pressure-volume interaction. 

Artificial production or dissipation by means of finite-difference 
approximations to the convective terms is still possible and must be 
eliminated for a valid simulation. kb shall again rewrite the convec- 
tive terms in the momentum equations in a different but equivalent form 
that not only prevents artificial kinetic energy production but regains 
the total energy conservation that we lost earlier. 

kb must ensure that the numerical method that we shall use is in- 
capable of artificially creating kinetic energy. To show this, we shall 
write the kinetic energy equation and integrate it over the periodic 
domain . 


Although we do not carry a kinetic energy equation in the simula- 
tion, its effect is implicit. The equation is formed from the mass and 
momentum equations . Therefore , what we do with these equations numeric- 
ally is reflected in the kinetic energy behavior. implying the chain 
rule to the time derivative of the kinetic energy shows 



UiUi 


where we see the appearance of 
equations . 


a “i'^i 3 

“i 3t ‘^“i " 2 3t 

l^pu^ and l^p 


p (2.2.1) 

for which we solve 


Integrating the kinetic energy equation over the domain, we find 
(with zero viscosity) 




It 


— = — dx + 


/ 3 ■»■ /* a 4- 

u. TT — pu.u . dx + / u. -» — P dx - 
i 3xj i j Jd i ^ 

-X 


Vi s 


( 2 . 2 . 2 ) 


yy — ^ — Pu . dx = 0 

2 3x . j 
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Because of the periodic domain, the integral of all quantities that 
appear totally inside a spatial derivative (divergence quantities) are 
zero. In addition, all terms that are evaluated on the boundaries, such 
as the constants of integration are also zero for the same reason. This 
is true both analytically and numerically, so those terms are not writ- 
ten here. 

Integration by parts shows that 

/ "i ^ ‘’"i X 4-^ slj ‘’V* ■ ° 

so that (2.2.2) becomes 

^ It ^ X "i 4^ 

Therefore, the only net contribution to the time development of the 
kinetic energy is the pressure term in (2.2.4). The convective terms 
have no net effect. Vfe need this behavior in the numerical simulation. 


Integrations are carried out numerically by sximmations. Mansour 
(1978) has shown that summation by parts holds for a wide variety of 
numerical methods. If we rewrite the convective terms as 


3 

Sx . 


PUiU. 


2 


•5T: 


+ U, 


3 

3x , 




(2.2.5) 


then summation by parts is valid and (2.2.3) is satisfied numerically. 
Vfe then, numerically, have the correct behavior for the kinetic energy 
that was described analytically in Eq . (2.2.4). 

As a side benefit, we have regained the conservation of total 
energy. Recall from (1.2.5) that the total energy is the sum of the 
Internal energy, (P/(Y“1)), and the kinetic energy, p(uj^u^/2). 
Therefore, the total energy behavior is determined by how we treat the 
mass, momentum, and pressure equations. Writing the time derivative of 
the total energy and integrating over the domain, we have 



3 ,•> 

3T P® 


i 3 


P dx + 



3 

It 


UiUi 

P-2— 


dx 


( 2 . 2 . 6 ) 
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Wfe have just discussed the numerical behavior of the last term. It 
Is represented by (2.2.4). The second term is obtained by integrating 
the pressure equation (2.1.10), resulting In 


X <-y- 


u §r + (r 


-1)^ 


^“l ■> 1 

- 

9x. 

1 


X“‘ 


9 

9x, 


P dx 


(2.2.7) 


Upon adding (2.1.10) with the appropriate viscous terms to (2.2.7), 
we find that the time derivative of the integrated total energy is 


r ’ 4. *"i H* 

- “i 157 "13 ^ 157 "ij 

Since numerical summation by parts holds, we find that 

P ^ = 0 




1 9x 


9X4 


C ^ 4 . 

I U , -5 T . . + -5 T . . dx 

i 9xj ij 9xj ij 


= 0 


( 2 . 2 . 8 ) 


(2.2.9) 


and therefore, numerically. 



3 

^pe dx 


0 


( 2 - 2 - 10 ) 


as it is analytically, so that we have regained total energy conserva- 
tion. 


This modification (2.2.5) to the equations is completely indepen- 
dent of the coordinate transformation. Vfe use both ideas in conjunction 
in the simulations. 
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2 .3 Summary of the Equations 

The previous two sections describing the coordinate transformation 
and the conservation properties provide the basis for the simulation. 
Since virtually all tensor equations that we shall discuss after this 
point have been transformed by the method of Section 2.1, we make the 
following definition. Recall from (2.1.5) that 


•HV ^ ®kj (2.3.1) 

Vfe shall absorb the into the derivative definition and drop the 
primes. After this point, when we write the transform is implicit 
in the derivative. For example, in the shear flow 
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3x^ 

^ - 


i = 2 




3 

3 x 3 

9 

tl 

This notation 
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be 

used unless 

specified 

otherwise . 

definition in 

mind , 

we 

present the 

equations as 

simulated 


(2.3.2) 


Keeping this 


Mass Equation 


Tt 


p + (pu^) = 0 


(2.3.3) 


Momentum Equations 


+ - -«iiS(PU2) + u[uj + 


(2.3.4) 


Pressure Equation 


^ p + S(U2_1 + 
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(2.3.5) 



where 


S = (mean velocity gradient) 

dy 

'ij ■ 

“jj ■ "'".JJ 

It should be emphasized that numerical considerations influenced 
the form that these equations (2.3.3) to (2.3.7) take, but they are in 
analytic form, with no approximations. The numerical method used to 
advance the flow fields in time is discussed in the next section. 

\ 

2.4 Desire for High Time and Spatial Accuracy 

These simulations are time-dependent and three-dimensional . It is 
necessary to choose ntunerlcal methods to compute the spatial derivatives 
and to follow the development in time. As stated in Section 1.5, we 
desire as accurate a simulation as possible, if the flow fields are to 
be used as a data base . 


(2 .3.6) 
(2.3.7) 


Before discussing the numerical methods, we must discuss the con- 
cept of stiffness . The one-dimensional Euler or inviscid Navler-Stokes 
equations serve to illustrate this property. The equations are 


° (2.4.1) 

^ (2.4.2) 

Tt ^ ° (2.4.3) 

Eqs. (2. 4. l)-(2. 4. 3) can be rewritten in vector form 

+ = 0 (2.4.4) 


where 
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(2.4.5) 
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u 

p 
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-1 

Q = 

u 

5 £ = 

u 
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p_ 


_0 

YP 

u 


A similarity transformation exists (Warming, 1975 ) that converts the 
Jacobian matrix, ^ to a diagonal form. Let ; 

P’^JP = A (2.4.6) 

where A is a diagonal matrix containing the eigenvalues of ^ and 
P is the corresponding eigenvector matrix. It is found that : 


A 


u 0 0 

0 u+c 0 

0 0 u-c 


YP 2 

— = c 
P 


(2.4.7) 


where u is the convective speed and c the speed of sound. This 
diagonalizjition is used on the nonlinear equation (2.4.4) only to show 
the various speeds at which information propagates . With a nonlinear 
equation it is not useful for the solution process, because the eigen- 
vector matrix P varies with time. 

There are several characteristic speeds at which information propa- 
gates : the convective and acoustic speeds and their sum and difference. 
In homogeneous turbulent flows, the fluctuating velocity, i.e., the 
convective speed u, is usually much less than the acoustic speed c. 
It is believed that the flow develops at a rate corresponding to the 
convective speed. However, in explicit numerical methods the size of 
the allowable time step is limited by events that develop at a rate 
corresponding to the acoustic speed. Hence, with these methods, numer- 
ical stability requires a very small time step (and consequently many 
time steps) so the cost will be high. 

This can be illustrated by solving the linearized form of (2.4.5) 
where is made constant by linearizing about the initial condi- 
tions. Ife subject the solution to periodic boundary conditions and a 
given initial condition Qq . 

Assuming that the equations have periodic solutions , we expand the 
solution in spatial Fourier transforms. Since this is a linear problem. 


20 



all Fourier inodes behave Independently. Therefore , for simplicity we 

choose a single representative mode associated with wave number k. 


Fourier transforming the linearized version of (2.5,5) results in 
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(2.4.9) 


which is a coupled set of linear ordinary differential equations with 
constant coefficients. 


life solve the system (2-5-9) by diagonalizing it in the fashion of 
(2-4-6). 

A AAA 

P“^ -Ir Q + P“^ J P p“^ ik Q = I- R + ik A R (2.4 .10) 

— o9t^ oooo ^ at o ^ ' 


resulting in an uncoupled set of ordinary differential equations . The 
exponential solutions of (2-5-10) are subject to the initial conditions 

A A 

R = P“^ Q (2 .4.11) 

o o o 

and the solution is transformed back by 


resulting in 


Q 


= P R 
o 


(2.4.12) 
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p — + — cos (c kt) - i — sin(c kt) J 
o '■ Y Y o - c o ' 

o 

Uq 

c ( — cos c kt - 1 — sin c kt) J 
o ^c o Y o ^ 

o 

Yu^ 

P^ [cos(c^kt) - i — ^ sln(c^kt) J 


(2.4.13) 


It can be seen that information propagates at different speeds in- 
volving Uq and Cq . If a simulation is to be accurate , it must 
properly represent all of those propagation speeds . The stiffness of 
the systems (2.4.4) and (2.4.9) are characterized by the wide range in 
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magnitude amongst the eigenvalues of A requiring very small time steps 
for ntimerical accuracy. This can be a severe limitation in that the 
speed of sound c^ can be many times the convective speed u^ . To see 
a significant development in the flow could take a tremendous amount of 
processor time. 

2-5 . Numerical Method Characterization 

kfe have spoken of explicit methods and their restrictive stability 

C 

criteria in compressible flows . In recent years implicit methods have 
been developed ( Whrming & Beam, 1978; Briley & ^facDonald, 1973) into 
very sophisticated algorithms that can greatly exceed the restrictions 
on time step required by expicit methods. 

Naturally, we are interested in such methods, as they could sig- 
nificantly reduce the computer time required. However, some procedure 
for comparing the behavior of the various methods has to be chosen. We 
choose the von Neumann analysis, (Lomax, 1967). This method is illus- 
trated by application to two very simple numerical methods , one expli- 
cit and one implicit. 

We choose a linear test equation, the linear form of the viscous 
Burger's equation. 

|_u+a|f= v^u (2.5.1) 

dX 

Wfe again Impose periodic boundary conditions and Fourier transfom the 
equation in space : 


^ u + (ika + vk^) u = 0 


The solution is : 

-(ika+vk^)t 

u = u e 
o 


At t = nh, the solution is; 


-“n 
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— (ika+vk )h 
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n 


(2.5.2) 


(2.5.3) 


(2.5.4) 
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Equation (2.5.4) shows that the solution at time (n+1) is related 
to time (n) by the constant factor: 


a 


-(ika+vk‘^)h _ -X 


= e 


12 13 

= ...(2.5.5) 


Ihe factor a holds the information about how the solution devel- 
ops in time and space. All numerical methods produce an approximation 
to this factor, and their performance can be evaluated and compared by 
how well they approximate this exact linear solution. 


The spatial derivative information is represented by the wave num- 
ber, k. All difference methods produce an approximation to k, so we 
may Include their effect on the numerical solution by introducing the 
"modified wave number" k' and using it to replace the exact k in the 
following analysis . fe show how k' is related to k as follows . 


For periodic solutions , the spatial derivatives can be represented 
in a general way by a discrete Fourier transform. Let 
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Take the derivative of (2.5.6) and equate it to, for example, a central 
finite difference. 
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Comparison of (2.5.7) and (2.5.8) shows them to be identical if 


k 

n 


sin k Ax 
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(2.5.9) 


Hence the concept of a modified wave number as used by Mansour (1978). 
The modified wave numbers for other difference schemes can be derived in 
a similar manner. Therefore, the characteristics of each type of dif- 
ferencing scheme can be analyzed in a general sense by using an equation 
like (2.5.10) to represent a derivative. 

ikx. 

f(xj) = ^ ik;^ f(k^) e (2.5.10) 
n=(-N/2) 

kfe can illustrate the behavior of some finite-difference schemes by 
plotting the modified wave number, k' , against the analytic wave num- 
ber, k. This is done in Fig. 2-1, where we show the behavior for 
second- and fourth-order central approximations to the first derivative . 

also show the pseudo-spectral derivative behavior where discrete 
Fourier transforms are used directly, to calculate the derivative. 

A perfect numerical derivative is represented by a 45® line on Fig . 
2.1. It can be seen that the behavior of the fourth-order method is 
better thjin that of the second order method, but both methods do not 
represent the high wave numbers or the small-scale structure well . The 
pseudo-spectral method is exact up to the maximum wave number represent- 
able on the mesh . Higher wave numbers in the solution are misrepresen- 
ted as lower wave numbers. They masquerade as contributions to resolv- 
able wave numbers; hence they are called aliased. With the pseudo- 
spectral method, care must be taken not to allow information into wave 
numbers greater than this maximum, or the solution accuracy will deter- 
iorate and instability can result. The pseudospectral method is far 
more accurate than the finite-difference methods, but is strictly lim- 
ited to periodic flows . These derivatives are the most suitable for our 
simulations . 
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Analysis for a Simple Explicit Method ; Euler 


The Euler explicit method represents the time advance. 

= u'^+h|^u” (2.5.11) 

Using Eqs. (2.5.1), (2.5.2), and (2.5.11), 

^nfl ''n ''n , ,,, 2. ^n 

u = = u - h(ika + vk ) u 

Ogg = 1 - (ika + vk^)h = 1 - X (2.5.12) 

It can be seen that (2.5.12) approximates (2.5.5) only to first order 
in X — hence Its classification as a first-order method. Other mea- 
sures of the method's behavior occur in dissipation error and phase 
error, as illustrated below. 

Let 


C = k ah ; V = (2.5.13) 

C and V are the Fourier space equivalents of, respectively, the Cou- 
rant number and the viscous stability number; they depend on both the 
time step h and the modified wave number k' . 

In Fig. 2.2, V has been set equal to zero. The two curves show 
the behavior of a and as the Courant number C increases. The 

analytic solution a neither decays nor grows in magnitude, but the 
Euler solution grows monotonlcally with time. It is unconditionally 
unstable for V = 0. The argument of Ugg can be at most ir/2, but 

the argument of a increases indefinitely. Phase information is 

increasingly inaccurate as C increases . 

Figure 2.3 illustrates the behavior of a and Ogg when C = 0. 
In this case the solution should simply decay. It can be seen that 
Ugg follows a fairly closely up to about V = 0 .3 and then diverges , 
and its magnitude finally becomes greater than 1 at V = 2.0, after 
which the method is unstable. 



Dissipation error (actually, antidissipation or growth) is evident 
in Fig. 2.2 in that the magnitude of Ogg does not remain 1.0, and in 
Fig. 2.3 in that the viscous dissipation is inaccurately represented for 
higher viscous stability numbers. For general X = (iC + V), these 
effects work in combination. 


Simple Implicit Method; Crank-Nicholson 

Wfe neixt present a simple implicit method. Certain advantages of 
implicit methods are shown by the Crank-Nicholson method. It will also 
be shown why methods of this type are insufficient for accurate time- 
dependent simulations when used with time step greater than the explicit 
stability condition. 


The (iirank-Nicholson method and implicit methods in general use the 
value of the unknown time derivative at the new time step. This in- 
volves solution of coupled equations which must be simultaneously 
solved, hence the designation "implicit." Crank-Nicholson is represen- 
ted by : 
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(2 .5.14) 


Analysis similar to (2.5.4) and (2.5.12) produces 


CN 


- 

(1 + j X) 


= 1 - X + -i- X^ - i X^ + ... (2.5.15) 


which approximates (2.5.4) to second order in X. Figure 2.4 shows the 
phase angle 6 of when V = 0 and Pig. 2.5 the magnitude of 

when C = 0. 

With V - 0, the magnitudes of a and are identically 

equal to 1.0, for all values of C. In other words, Crank-Nicholson 
produces neutrally stable solutions for any size time step. Fig. 2.4 
plots the arguments of a and “gjg* It can be seen that, while Crank- 
Nicholson solutions are stable for Courant nimbers greater than 1.0, 
their accuracy as represented by phase infomation can be very poor 
(likewise for other implicit methods). 


26 



Figure 2 .5 shows that for C = 0 the solution can be alternately 
positive and negative at large viscous stability numbers because Uqjj 
becomes negative. While this behavior Is stable, It Is an Incorrect 
representation of the true behavior of the solution (2.5.4). 

Vfe can apply this type of reasoning to more complicated sets of 
test equations such as (2.4.9) that have several propagation speeds. 
One finds that, If all motions In the problem are to be simulated accu- 
rately, the Courant number based on the fastest propagation speed must 
be of the order of magnitude of 1.0 or less. Also, the viscous stabil- 
ity number must be 1 .0 or less . 

Since implicit methods require more work per time step, their re- 
striction to small Courant and viscous stability numbers makes them 
inefficient . Vfe therefore abandon implicit methods for our time-accurate 
simulations . 

Some work was performed in incorporating the previously discussed 
conservation properties into an implicit method. Additional problems 
were encountered in designing an efficient implicit algorithm that in- 
corporated highly accurate spatial derivatives. As stated above, the 
implicit methods are not suitable for highly accurate time-dependent 
simulations. They lack sufficient accuracy for a time-accurate solu- 
tion. They may find use in simulations of variable density, low Mach 
number flow, where the acoustic speeds are unimportant in the solution, 
but are troublesome numerically. 

2 .6 Advanced Explicit lyfethods 

There are two fully explicit methods that were finally considered. 
MacCormack's (1969) method was analyzed and programmed into a three- 
dimensional code to run on a CDC 7600. The final method of choice, 
however, was a combination of the classical fourth-order Runge-Kutta 
with spatial derivatives based on fast Fourier transforms, the pseudo- 
spectral method. 

MacCormack's method is well known, having been successfully used in 
many compressible flow codes over the last decade. It is a modification 
of a Lax-Wfendroff scheme that incorporates forward and backward differ- 
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ences alternately on the predictor and corrector steps. Its formal 
accuracy is second order in both space and time. The reader is referred 
to the original paper , in which the method is laid out and estimates 
made of its accuracy and stability. 

For purposes of comparison, we shall apply the analysis developed 
in the previous section to compare these two methods . 

In most numerical algorithms for partial differential equations, 
the method for representing the spatial derivatives can be discussed 
separately from the method for the time advance . To a limited extent , 
one is free to choose combinations. MacCormack' s method, however, is a 
specific combination of the two and derives its unique qualities from 
this combination of methods . 

Following the analysis in Section 2.5, the computational root of 
MacCormack's method is 

a = 1 - V - y C'^ + i V'^ + i(C" V’ - C” ) (2.6.1) 


where 


V = hvk'^ ; C = hak’ ; C" = hak" (2.6.2) 
c 


and 


k' 


sin k 


Ax 
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L-x-H 


k” 


sin kAx 
Ax 


(2.6.3) 


The forward and backward differences result in the two modified 
wave numbers of (2.6.3). These wave numbers correspond to central 
difference approximations to the first derivative that are taken over 
mesh spacings of Ax/2 and Ax, respectively. The highly dissipative 
behavior of this method results from the combination of modified wave 
numbers (2.6.3). 

Figure 2.6 shows the computational root for zero viscosity (V = 
0). Both time and spatial behavior are represented. Keeping in mind 
that the exact solution lies on the unit circle, it can be seen that for 
high wave numbers (kA approaching x) or large time steps, the 
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dissipation is severe without the inclusion of viscous terms . If an 
initial flow field containing all resolvable length scales is advanced 
by this method, the smallest scales (highest wave numbers) will swiftly 
disappear through the action of the numerical method. 

It was found in an isotropic homogeneous flow that this numerical 
viscosity was sufficient to maintain a stable calculation at a Taylor 
microscale Reynolds number of several thousand . At the relatively small 
mesh size of 16x16x16 = 16 , this is unphysical* In addition, 
statistics that have to do with high wavenumbers , such as velocity de- 
rivative skewness , were incorrect in comparison with experiment . 

Since an accurate method for all resolvable wave numbers is 
desired, this method must be abandoned for our purposes. 


2 .7 Method of Choice 

The numerical method that we have chosen has a combination of the 
highest time accuracy that can conveniently be Implemented on the ILLIAC 
IV and the highest spatial accuracy. Ife shall use fourth-order Runge- 
Kutta for the time advance and the pseudo-spectral method for the spa- 
tial derivatives. For completeness we show both. 

Fourth-order Runge-Kutta is written: 
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and the pseudo-spectral spatial derivatives directly implement 
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n n 


(2.7.2) 


n=(-N/2) 

by means of the Fast Fourier Transform algorithm. 
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The computational root for pseudo-spectral fourth-order Runge ^tta 
is, by the kind of analysis used above. 

1 \ j. 1 i ,3 , 1 ,4 /o 7 o\ 
ot * i”X+'2*x ^ 

For zero viscosity (v = 0), this root is depicted in Fig. 2.7. 
For Courant numbers less than one, the method is very accurate, since it 
closely follows the analytic root. In using this method, the time step 
is chosen such that the Courant number associated with the maximum wave 
number is one. Courant numbers associated with the lower wave numbers 
are then less than one, and the computational root more accurately fol- 
lows the analytic root for these larger scale motions . 

Figure 2.8 shows behavior of the viscous terms alone (C = 0). In 
this combination of methods, the computational root accurately approxi- 
mates the analytic root up to a viscous stability number of one. 

The combination of fourth-order Runge-Kutta and pseudo-spectral 
spatial derivatives requires a large number of operations (and conse- 
quently CPU time) per time step. Additionally, it requires four levels 
of storage in the computer for intermediate predicted steps, placing a 
burden on the memory resources of the machine for large meshes. This 
compares to two levels of storage for second-order methods. However, we 
are able to use the method at a Courant number of one, in contrast to 
many common second-order explicit methods that must have Courant numbers 
much less than one. Although the method requires more CPU time per time 
step, its time steps are correspondingly longer than any suitable 
second-order method . 

We are still limited by accuracy constraints to a Courant number of 
one. To perform a complete simulation in a compressible flow with the 
large mesh system that we wish to use will absorb an incredible amount 
of processor time. In the next chapter we show how this method is used 
to solve Eqs. (2 .3 .3)-(2 .3 .7) in a code for the ILLIAC-IV computer at 
NASA- Ames Research Center. Wfe shall discuss how the code is constructed 
and tested and show the amount of CPU time that it absorbs. 
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Chapter III 


PROGRAM DEVELOPMENT AND TESTING ON THE ILLIAC IV 


The construction and testing of the computer code are described in 
this chapter, as well as the characterization of a typical run. 

3.1 Ihe ILLIAC IV 

The ILLIAC IV is a large, one-of-a-kind, parallel, processing 
computer built by Burroughs. Codes designed for the ILLIAC IV can be 
run only on it (although the algorithms could be converted to other 
vector computers) . Illiac IV is composed of 64 individual computers 
called PEs that operate in lockstep under the control of a central 
managing computer called the CU. The PEs have access to the large 
rotating disk memory that contains the bulk of the memory on the 
computer. The in-core memory of each PE is quite small, so all of the 
flow fl€ild data resides on the disk and is brought into core a small 
piece at a time to be operated upon. The efficiency with which this 
transfer is accomplished has a great effect on how fast a code will 
run. This operation is directed by a set of Instructions called a disk 
laap that controls the transfer of specific data between the disk memory 
and core. 

3.2 The Pencil Data Management System 

The flow field must be divided into regions that are transferred 
into core one at a time. Because the fast Fourier transform is used so 
extensively in this code, data from a line of mesh points that extend 
entirely across the computational box must be in core all at once, the 
memory-management system chosen is the "Pencil System," as developed by 
Pulliam and Lomax (1979). 

The ILLIAC IV is run in 32-bit word mode, which allows the pencil 
size to be 8 x 16 x (mesh size) words in size. This is depicted in 
Fig. 3.1. The spatial derivatives are handled, in order, such that 
all derivatives in the x-direction are performed while the pencils 
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in that direction are in core. The other two directions are handled 
consecutively, in a similar manner. A typical predictor or corrector 
step involves processing each of the x pencils in turn, then each of 
the y pencils , and then the z pencils . Each set of pencils requires 
one complete pass over the entire flow field . Since fourth-order Runge- 
Kutta requires four function evaluations per time step, the code 
requires twelve passes over the flow field . While this is very time- 
consuming, the pencil system allows great freedom in constructing means 
of analyzing the simulated flow fields and allows easy modification of 
the code for other purposes. 


3 .3 Coordinate Remeshing 

The linear coordinate transformation as described in Chapter II is 
analytic and is incorporated into the equations in its entirety. This 
transformation allows the transformed coordinates x and y to become 
almost parallel after a length of time. Tb avoid accuracy problems in 
the code, it is necessary to stop the calculation after a time, shift 
the sheared mesh back to the other side of Cartesian, and then proceed 
again (Fig. 3.2). It has been found that this remeshing process is 
essentially exact if the calculation is stopped when the top of the 
computational box has moved one-half period (St = 0.5). The mesh is 
then shifted back one period to (St = -0.5). When the remeshlng is 
done at this point, there is no Interpolation necessary, as the points 
of the new and old meshes fall exactly on top of each other. There is 
one numerical problem associated with this . Because the orientation of 
the mesh in the y-directlon is changed, there is aliasing in that 
direction associated with remeshlng. Let a particular Fourier mode be 
represented in two dimensions on the second mesh in Fig. 3.2 as 


u(x,y) 


u(kj^,k^) 


i(kj^3cfk2y) 


(3.3.1) 


(Recall that the mesh system shears with the mean velocity .) If the 
remeshing is done at St = 0.5, then the new coordinates, as measured 
on the third mesh in Fig. 3.2, are 
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X + y X 


(3.3.2) 


r* - 


y» = y 


Substituting (3-3-2) into (3-3-1), we find that the same Fourier mode is 
represented on the "remeshed" mesh as 


A ik, x' i(k,-k,)y' 

u(x'.y') = u(kj^,k 2 ) e e 


(3.3.3) 


If (k 2 ”k 2 ^) falls outside the range of available wave numbers, that 
information will appear in a lower wave number in the y-direction. This 
problem is eliminated by filtering the high wave number parts of the 
solution in the y-dlrection both before and after the remesh process . 
The filtering is performed by Fourier transforming all the flow vari- 
ables in y, truncating the top one-third of the wave number coeffici- 
ents and retransforming back to real space . This is performed on flow 
fields which contain little information in the highest wave numbers and 
only at remeshings . It removes at most 1% of the turbulent kinetic 
energy. Its effect on the solution is small, and it eliminates a known 
source of error . 


3 .4 Initial Conditions 

Turbulent initial conditions as used in low Reynolds number and 
large eddy simultions will always be somewhat artificial. There are 
many statistics used to describe these flows but no organized algorithm 
for producing these statistics in an initial flow field. The hope, that 
is generally borne out by numerical experiment, is that, if some statis- 
tics are enforced on the initial field, the others will develop through 
the action of the equations. In other words, the initial flow fields do 
not represent a turbulent flow field, but after being advanced for a 
time by the code they develop the characteristics that allow us to call 
them truly turbulent . 

In the compressible flow, the five flow variables are completely 
independent and subject only to the restriction that both density and 
pressure be positive. 
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The variables may be thought of as being In nondlmenslonal form . 
They are non-dlmenslonallzed on a sound speed c^, a density Pq> 
and the computational box length . Because of this nOndlmenslonal- 
Izatlon the velocities are closely related In magnitude to Mach num- 
bers. The average value of the computational density Is one, as Is the 
computational pressure. This makes the average Initial computational 
sound speed c (different from c^) of magnitude 1.18. 

c^ = Y^, c = /y = 1.18 (3.4.1) 

Ife have chosen to start each calculation with a velocity divergence 
free field, as In an Incompressible simulation. Since the Initial flow 
fields are artificial, any dilatation Introduced Is unphyslcal . If any 
Is to appear. It should grow through the action of the equations of 
motion. Wray (1980) has found In two-dimensional simulations of Iso- 
tropic turbulence that these Initial divergences decay very quickly, 
leaving behind density variations and a flow that behaves almost Incom- 
presslbly. In the shear flow we felt It best to allow compressible char- 
acteristics to develop under the forcing action of the mean shear rather 
than Introducing them as Wtay did in the unforced isotropic flow. 

The procedure for producing the Initial velocity field is as fol- 
lows . Ivb first choose a completely arbitrary random set of velocities 
through the action of a random number generator. Vfe add the gradient of 
a yet to be determined potential function onto this velocity field. 

''l " (3.4.2) 

where u® is the solenoldal (dilatation-free) velocity field and u^ 
Is the original random field. Since u| is solenoldal, it disappears 
when we take the divergence of (3.4.2), producing a Poisson equation for 



Because of the periodic boundary conditions and the availability of 
fast Fourier transforms, we have a particularly convenient and accurate 
way to solve this equation. It is^ first Fourier transformed. 

9 ^ 2 

- k <|) = - ik.u, (k = k.k ) 

j J 3 J 

" ik.u 

<!> = (3.4.4) 

k 

Inverse Fourier transforming provides the solution (j>. The con- 
stant of integration is fixed by setting the zero wave number coeffi- 
cient to zero . 

The velocity field u| is solenoidal, real, and free of mean vel- 
ocity but has a white noise type of velocity spectrum. fe take this 
spectrum by transforming the velocities into Fourier space and integrat- 
ing the kinetic energy in spherical shells to produce a 3-D energy spec- 
trum . 


E(k) = < u^(k) u*(^) > (3.4.5) 

Here, < > indicates an integration over all directions, leaving the 
kinetic energy as a function of the wave number magnitude alone. 

All of the velocities in each spherical shell are then adjusted by 
the same constant (a function of k) to enforce a specified energy 
spectrum, E(k) onto the flow field. 

Wb silso enforce that the time derivative of the divergence be zero, 
as in an incompressible flow, so that there is no violent behavior when 
the simulation starts. This is done by specifying the pressure field 
through another Poisson equation. take the divergence of the momen- 
tum equations 


83 .3 8 . „2 

pu, + - 5 ;^ pu. u, + V p 


3x. 8t ^ i 3x. 3x, *^“i“j 


■5^ "Sir ij 


(3.4.6) 


For the initial flow field with zero dilitatlon and constant density, 
(3.4.6) becomes 
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(3.4.7) 


V^P = 


which we solve in the same manner as 


9u^ 9uj 
J i 
(3.4.3). 


All flow fields used as Initial conditions were set up In t|ils man- 
ner. Initially they are incompressible, but their development is then 
guided by the full Navier-Stokes equations. Velocity divergences de- 
velop quickly. 


3 .5 Performance of the Code 

There are actually three codes involved in the simulations . Ihe 
severe limitation on in-core memory of each processing element (PE) has 
forced the dissection into codes : 

1. TBr-LOD.IHICS . Creates and stores the initial flow fields. 

2. TB-LOD.HSTA. Ihe time advance code, which does some prelimi- 
nary data reduction. 

3. TB-LOD.REDUCE-DATA. Performs the bulk of the data reduction. 
The mnemonics TB is the author’s identifier on the ILLIAC system. The 
LCD indicates an executable load module. IHICS stands for isotropic 
homogeneous initial conditions. HSTA means homogeneous shear time ad- 
vance . 

All three codes operate with the same pencil data-management sys- 

3 3 3 

tern. This system can operate with three mesh sizes, 16 , 32 , or 64 • 
For the time-advance load module, the code requires 28 words of memory 
for each mesh point. In 64 , this is 7.34 million words. Including 
additional disk areas needed for data output and temporary storage, the 
management system allocates 12.17 million words out of an available 15.9 
million. Although the data-reductlon load module does not require the 
large amount of memory that the time-advance module needs, this storage 
capability is very convenient in processing the simulated flow fields. 

The time-advance load module consumes most of the central processor 
time. In 32^ form, this load module requires 8.5 seconds per time 
step. In 64^ form, it requires 89 seconds per time step. Most of the 
work was performed with the 64^ version. 
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There are eight simulations associated with this thesis . The 
shortest, corresponding to the highest shear rate, took 7.5 hours. The 
longest, the isotropic simulation, took 18.0 hours. This does not in- 
clude all the data reduction, which was also performed on the Illiac. 


3.6 Testing the Code 

As was stated in Oiapter II, one must be particularly careful to 
run checks on the code to ensure that it is properly constructed. To 
this end, checks were performed in three categories. The code was 
written to be quite readable and was thoroughly sight checked. The 
second category involved simulation of two-dimensional Taylor-Green 
problems, which are known to be exact solutions to the incompressible 
flow equations . The third category involved simulation at low Mach 
number of isotropic turbulence and comparison of some statistics with 
experiment . 


3.6.1 TVo-Dimenslonal Thylor-Green Problems 

There is an exact analytic solution to the two-dimensional, incom- 
pressible, Navier-Stokes equations. It is known as the two-dimensional 
Taylor-Green solution. The solution is as follows: 

1/2 -(kj^+k^)^! 

u = - A cos kj^x sin k^y e 


V 


a ^/2 


-(kl+kjjvt 

Sin lex cos k2y e 


(3.6.1) 


A -2(kJ+k^)»t 

P = ~ 7- (cos 2k, X + cos 2k.y) e 

4 1 2 

The flow fields consist of 2-D vortices, arranged rectangularly, that 
simply decay in strength with time. 

This is not an exact solution for the full compressible equations. 
It was thought that, for low MEich numbers, the compressible flow field 
should behave almost Incompressibly and the solution (3.6.1) should be 
closely approximated by the compressible code . By numerical experiment 
this is true . 
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The initial conditions were Eqn. (3.6.1) at t = 0, “ ^2 ~ 

and A = .05 (about Mach = 0.04). The code was first run with a Car- 
tesian mesh (St = 0) . After 40 time steps the maximum deviation from 
the analytic solution was 1 .42 x 10~^ or .028% of the maximum veloc- 
ity. KLth the mesh sheared to St = 1.0, which preserves the periodic 
boundary conditions on the solution, the corresponding numbers are 
1.47 X 10"^ and 2.94%. 

The numerical solution is in quite close agreement with (3.6.1), 
even though it is not an exact solution of (1 .2 .2)-(l .2 .7) . These 

o 

experiments were run in 16 and for all orientations of the initial 
conditions in x, y, z for the St = 0 mesh. For St = 1.0, they were 
performed in the x,y plane (downstream, cross-stream) . 

3.6.2 Low Mach Niimber Isotropic Turbulence 

As in the previous section, it is felt that at low Mach numbers the 
flow fields should be virtually incompressible and therefore should pro- 
duce turbulence statistics quite similar to incompressible codes and to 
experiment. Again this is true by numerical experiment. 

O 

A 64 run was performed with the shear rate set equal to zero. 
The initial conditions were constructed by the code TB-LOD.IHICS in the 
manner of Section 3.4. The initial spectrum was a box between wave num- 
bers 8 and 16. The initial average Mach number was M^ - .078 and the 
initial Ihylor microscale Reynolds number was Re^ = - ^0* 

Figures 3.3 to 3.10 show the evolution of the flow^cluring the simu- 
lation. Figures 3.3 and 3.4 are three-dimensional energy spectra. In 
each figure the lower curves are 3-D spectra of the normal stresses of 
the Reynolds strss tensor. It can be seen that the flow evolves from 
the very artificial initial spectrum to a realistic-looking low Reynolds 
number spectrum as the simulation proceeds . There is no linear region 
to represent an inertial range, but the slope does pass through the 
value of -5/3 as the wave number Increases. At low k the slope 
passes through a value of 4, an analytic shape for low wave numbers. 
After wave number 32, a steeper slope in the spectrum is observed. This 
is a result of the way the spectra are taken. The scalar wave number 
k is the magnitude of the vector k. The energy in each spherical 
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shell Is simply stunmed over the number of points (samples) in that 
shell. After wave number 32, the shell begins to grow outside the cubic 
domain of k and therefore the number of contributing points decreases 
above k. = 32. The decrease in energy simply reflects the decrease in 
sample size . No attempt is made to extrapolate beyond k = 32 . Figure 
3.5 is representative of lateral and longitudinal two-point correlations 
in the ai-direction; z was arbitrarily chosen for display, as all di- 
rections are equivalent . The figure was produced at the same time as 
Fig. 3.4, t = 7.8. It shows the expected primarily positive behavior 
of the longitudinal correlation R 33 ( 0 , 0 ,r 3 )» 




< u^(x) Uj(^r) > 

< u^(x) Uj(x) > 


(3.6.2) 


and the close agreement between the two lateral correlations , 


R, 


22 ( 0 , Ojt^) and Rj^j^(0 ,0,r2) • Figure 3.6 shows the three longitudinal, 

Ej^j^(kj^), E22(k2)» ^33^^3^’ 


one-dimensional spectra, 
defined as 


E^^(k) = u*(S) dk 2 dk^ 


(3.6.3) 


®22^^2^ and E 23 (k 2 ) are defined similarly. The similarity of the 
three curves shows that the flow remains quite isotropic at all wave 
numbers, as it should. The one-dimensional spectra are Fourier trans- 
forms of the two-point correlations. Because there is little informa- 
tion in the high wave number region, microscales as calculated from the 
correlations should be quite accurate. Figures 3.7 and 3.8 show the 
time evolution of the microscales and the Integral scales , respectively . 
The Taylor microscales are defined in the usual way by fitting the 
osculating parabola to the two-point correlation at small separation. 
The integral scales are defined as twice the separation where the two- 
point correlation first reaches a value of 0.1. The reason for this 
will be explained in the next section. After t = 5.0, both micros and 
integral scales grow almost linearly, but very slowly. 

Figure 3.9 shows the time evolution of the kinetic energy and its 

three nomal components . A power-law fit to this curve shows a best 

fit with a slope of about -1.25, in good agreement with theory (Hinze 
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(1975)) and experiment (CBC (1971)) for low Reynolds number isotropic 
turbulence . 

Figure 3.10 shows the time development of the velocity derivative 
skewness defined as : 

3u. 3 / 3u. 2 3/2 

Ski ~ ^ ^ ^ summation) (3,6,4) 

The skewness is associated with the energy cascade process and is mea- 
surable in the wind tunnel for the downstream direction. Its value is 
quoted in several experiments (Tavoularis et al. (1978)) as about -0.40. 

It can be seen that all three curves in Fig. 3.3.0 come to as3mip- 
totlc values of about -0.35 to -0.41. 

Although the time-advance code is designed for compressible flow, 
these results give confidence that it is capable of simulating low Mach 
number (Incompressible) turbulence. It must be able to simulate this 
flow before it can be used in a higher Mach number problem. 

3 .7 Description of a Typical Sheared Run 

A typical sheared run is described, along with the limitations and 
troubles generally encountered. This is presented to show the limits of 
validity of these simulations and the criteria used to judge this valid- 
ity. 

As described in Section 3.4, the initial conditions used to start 
the simulation are quite artificial. The initial conditions are simply 
described as a constant density field, a dilatation-free, random, iso- 
tropic, initial velocity field with a square wave (top hat), 3-D energy 
spectrum, and a pressure field set according to the velocity field to 
maintain 3/3t( 3u^/3x^) = 0 at time zero. As described in the previous 
section, the flow field requires a certain time to develop turbulent 
flow field characteristics. In the isotropic flow this was judged as 
the time when the velocity derivative skewness reached an asymptotic 
value . 

In the shear flow, this time was judged by several different stan- 
dards. Fortunately, they all point to a common time when the flow field 
might truly be judged turbulent . These criteria are ; 
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1 . 


The shear stress coefficient reaches an asymptotic value. 

< u'v' > /V< > < v'^ > = const (3.7.1) 

2 . The kinetic energy is seen to Increase . 

3. The 3-D energy spectrum reaches an asymptotic "shape" as it 
changes from the initial square wave spectrum. 

Ihe last criterion is a matter of judgment, but estimates made by cri- 
terion 3 agree quite well those with criteria 1 and 2. 

After this initial period, during which the flow field evolves from 
the artificial initial conditions, there is a period during the simula- 
tion when the judging criteria indicate that a true turbulent flow field 
is being simulated. These criteria will be discussed below in the de- 
scription of the simulation HS64B. This period, during which the sta- 
tistics are extracted from the flow field, is thoroughly discussed as 
the subject of Chapter IV. 

The simulation gradually moves out of this valid stage into a phase 
in which the scales of motion grow too large for the computational box. 
In the homogeneous shear flow, the length scales associated with the 
turbulent motion are known to grow with time. This indicates that 

structures or eddies in the flow grow to a size where they are influ- 
enced by the imposed periodic boundary conditions. Beyond this point, 
the flow field is no longer representative of an infinite, homogeneous 
shear flow, so the simulation is stopped. Tills time is judged by the 
appearance of the two-point correlations, which are a statistical mea- 
sure of the spatial relationships in the computational box. These 
correlations were defined in Eq . (3.6.2). 

As will be discussed in Chapter IV, the time is nondlmensionalized 
by multiplying it by the shear rate. In terms of this nondimensional 
time, St, all of the simulations were judged to be valid between 
times St =4.0 and St = 6.0. 

The simulation chosen to illustrate the history of a sheared run is 
labeled HS64B. As discussed in Section 3.4, the flow-field variables 
are nondlmensionalized on a density a sound speed Cq, and the 
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computational box length . Other nondlmenslonal quantities may be 
formed directly from the nondlmenslonal flow variables, in which case 
the original nondimensionalizing values of , and drop 

out. In these terms the initial characteristics may be summed up: 

1. Shear rate S = 1.5. 

2. Initial fluctuating velocity q = 0.1. 

3. Square wave 3-D spectrum 8 < k < 16. 

4. Box Reynolds number Re^j = C^Lp/v = 500. 

Some nondlmenslonal measures : 

1. Taylor microscales Reynolds number Rej^^ = 25.5. 

2. Average fluctuation Mach number M = .078. 

The initial spectra appear exactly the same as Fig. 3.3, so they are not 
repeated here. Figure 3.11 shows the development of the coefficients 
associated with the off-diagonal elements of the Reynolds stress 
tensor. The < u'v' > coefficient was defined in Eq . (3. 7.1). In this 
shear flow the other two coefficients associated with < u’w' > and 

< v’w' > should remain zero. It can be seen that this is so to a few 
percent through St = 6.0. The shear stress coefficient associated with 

< u'v' > starts from its isotropic value of zero at St = 0.0 and 
reaches an asymptotic value of -.64 by St = 4.0 and remains virtually 
flat through St = 6.0. 

The development of the kinetic energy as a function of St is de- 
picted in Fig. 3.12. It can be seen that it reaches a minimum at about 
St = 4.5 and Increases afterward. The < pu' > component of kinetic 
energy reaches a minimum much earlier at about St = 2.0, but energy is 
quickly drained from this term into < pv' > and < pw' > through 
the action of the pressure-strain terms in the dynamic equations for 
these quantities. 

The mechanism of kinetic energy production is obviously in opera- 
tion almost from the start of the calculation, but the shear stress co- 
efficient which represents this mechanism reaches its asymptotic value 
about the same time that the kinetic energy begins to Increase. 
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Figures 3.13 and 3.14 show the 3-D energy spectra at times St = 
4.0 and St = 6.0, respectively. It can be seen that they have 
evolved considerably from the initial spectrum and that they closely 
resemble each other. Kinetic energy originally present in the initial 
spectrum and the kinetic energy added to the flow through the action of 
the production term have been redistributed among the wave numbers in 
the 3-D spectra through the action of the nonlinear convective terms. 
Tlais process fills in both the low wave numbers, k < 8, and the high 
wave numbers k > 16. Ihe viscosity chosen for the simulation plays a 
great role in determining the shape of this spectrum (alternatively, the 
box Reynolds number) . The nonlinear convective terms will continue to 
propagate kinetic energy to higher and higher wave numbers (smaller and 
smaller scales) , unless there is sufficient viscous dissipation to 
change this kinetic energy into heat. If there is insufficient viscos- 
ity or, equivalently, the box Reynolds number is too high, the higher 
wave numbers collect too much kinetic energy, as it is not being 
dissipated quickly enough. The phenomenon of aliasing occurs, where 
information destined for higher, non-existent wave numbers returns to 
masquerade in, and pollute, the low wave numbers. This occurs when 
there is any information in the wave numbers greater than 2/3 of the 
maximum wave number. Figure 3.15 shows the one-dimensional, 
longitudinal kinetic energy spectra at time St = 6.0. It can be seen 
that information contained in each direction in the 2/3 wave 
number (in our case k = 21) is two orders of magnitude less than the 
energy peak at low wave number . 

The simulation begins to degrade when the structures in the flow 
grow large enough to be affected by the periodic boundary conditions. 
This is shown in Figs. 3.16 through 3.18. Figure 3.16 shows the two- 
point correlation at St = 4.0. Figure 3.16 shows the correlation at 
time St = 6.0 and Fig. 3.18 at St = 7.0. At St = 4.0, the curves 
still show almost zero correlation at Delta R = 3.2, which is half of 
the computaltonal box width. For a valid use of periodic boundary con- 
ditions, motions in regions separated by half the box width must be 
uncorrelated. By St = 7.0, the R33(0,0,r3) correlation has 
reached 18% at the half width, indicating Interference of the boundary 
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conditions . By limiting the maximum tolerable half-width correlation to 
10%, we restrict the region of validity to St < 6.0. 

In summary, each calculation passes through three stages: 

1 . Develoment away from the Initial conditions . 

2. TXirbulent simulation. 

3. Interference of the boundary conditions. 

The validity of the second stage is affected by the value of the box 
Reynolds number. There must be sufficient viscosity to hold the energy 
in the high wave numbers much below that in the peak wave number . For 
the simulations in this thesis, we use one order of magnitude. 

Ife have confidence that the computer codes are operating properly 
and now proceed to use them as a numerical wind tunnel . 

In the next chapter we discuss the seven complete compressible, 
homogeneous, shear flow simulations and present the measurements that 
have been made from the simulated flow fields . 
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Chapter IV 

CALCULATIONS AND RESULTS 


The calculations and results are discussed In this chapter. The 
independent dimensional parameters of the flow are listed, as are the 
nondlmenslonal groups formed from them. A description of the simula- 
tions that were performed is presented. Ivfe then discuss the behavior of 
the Reynolds stresses and the Reynolds stress equations, with particular 
attention paid to modeling of the various terms . 


4.1 Dimensional Parameters and Nondimensional Groups 


There are eight Independent dimensional parameters that can be used 
to characterize compressible, homogeneous shear flow. They are: 


1 . 

2 . 

3. 

4. 
5 

6 . 

7. 

8 . 


< u^u^ > 


dU/dy 


Density 

Moleculcir viscosity 
Turbulent intensity 
Speed of sound 

Mean velocity gradient (shear rate) 

Integral length scale used throughout) 

Taylor microscale 
Coefficient of heat conductivity 
Some of these quantities are related through the dynamics of the flow 
and are, therefore not truly independent. 


P 

P 

c 

S 

L 

X 

K 


The density is the average density in the domain, a constant. The 
molecular viscosity is fixed for a given simulation. The turbulent 
intensity is the trace of the Reynolds stress tensor divided by the 
density. The speed of sound is defined by = c . For a perfect 
gas, c “ p » where the ratio of specific heats Y is 1.4. The shear 
rate can be chosen arbitrarily. 

We define the integral length scale, L, as twice the separation 
at which the two-point correlation first reaches a value of one-tenth. 
This definition differs from the standard integral definition. The two- 
point correlations can have large negative loops, and therefore definite 
integrals of these functions can be poorly behaved. Figure 4.1 shows 
three correlations for a typical flow field. Harris et al. (1970) 
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encountered these problems and defined L by integrating up to the 
first zero crossing. His definition and ours are both portrayed. Our 
definition gives a length scale which is always larger than both Corr- 
sin's and the standard one. It has been adopted for ease of calcula- 
tion. The Taylor microscale is defined in the usual way by fitting the 
osculating parabola to the two-point correlation at small separation. 
The coefficient of heat conductivity is set by fixing the Prandtl 
number, using y = 1 with the perfect gas law constant. 

Dimensional analysis shows that five nondimensional groups may be 
formed. Three of these are shown in Table 4.1, where we also include 
values of these quantities found in our simulations and in experiments . 
Ivb choose values of the shear number that are similar to those in high- 
speed mixing layers and homogeneous flow experiments . Ub choose a range 
of fluctuating Mach numbers which spans the range from incompressible 
flow to values observed in high-speed flows . >fesh resolution limits the 
Taylor microscale . Ife use the highest value allowed . The ratio of 
length scales is known to be a function of the Reynolds number. It is, 
therefore, not independent. In our simulations it has a value of about 
4. As just stated, we have fixed the Prandtl number. In all our simu- 
lations we use Pr = 0.74, a value suitable for air. This leaves just 
three truly independent nondimensional groups that we vary during the 
simulations. 

Other nondimensional groups can be formed from combinations of the 
groups in Table 4.1. For instance, SL/q times q/c is the "shear 
Mach number," SL/c. Wb mention this number because it is more akin to 
the conventional external flow Mach number. It represents the change in 
Mach number across a typical large eddy. Note we take 

In some results we show the time-dependence . The nondimensional 
time is St . 

4.2 Description of the Simulations Performed 

Eight complete simulations were performed . Including the inevi- 
table waste, approximately 250 hours of central processor time were 
consumed on the ILLIAC-IV. Seven of these simulations were homogeneous 
shear flows. The eighth was the low Mach number isotropic flow that was 
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described in Chapter III. Wfe present data from the seven shear flow 
simulations in this chapter. As discussed in Chapter III, we start the 
code with artificial initial conditions, and the flow field must develop 
into a truly turbulent field through the action of the equations of 
motion. Consequently, in each simulation, there is a time before which 
we cannot accept the data. Correspondingly, there is a time after which 
the scales of motion are too large for the computational box and peri- 
odic boundary conditions are no longer appropriate. fe find the flow 
fields to be satisfactory between nondimensional times, St, of four 
and six in all runs, which we hereafter designate as the "good times." 

Turbulent statistics are evaluated only when St is an integer. 
Tills was done as a matter of convenience and because a short sample time 
would give flow fields that were not sufficiently independent. This 
leaves us with 21 flow fields to analyze. 

Figures 4.2 and 4.3 show the 21 flow fields in the parameter space 
of shear number, S, fluctuating Mach number, M, and Taylor micro- 
scale Reimolds number, 

In 'Dable 4.2 we tabulate the variation of the three nondimensional 
parameters and show our labeling for each of the simulations . there is 
no simulation labaled HS64E (there was, but it had far too high a vis- 
cosity and was discontinued) . In many cases we show plots for each of 
the seven simulations (seven plots), and we shall designate the plots by 
labels a through h. To make the plot numbers correspond to the flow 
labels , we always skip the label e . 

Initial data reduction was performed on the ILLIAC-IV. Calcula- 
tions of spectra and integrations over the 64 x 64 x 64 mesh require the 
entire flow field and are practical only on this machine. Hawever, once 
these qxiantities are calculated, they are transferred to a more 
conventional serial computer for further processing . 

The appendix presents the raw data from the ILLIAC and may be use- 
ful for further investigations . All correlation data that we present 
are further reduced from the appendix. 
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4 .3 Averaging for Statistics; Connection between Favre and Conventional 
Averages 

There are many types of averaging that have been used in turbulent 
flows . Vfe have presented the averaged equations and our approach to 
solving the unaveraged equations, but the averaging process is still 
unspecified. 

The ideal averaging process is an ensemble average, in which many 
realizations of the same flow are averaged. A homogeneous flow is sta- 
tistically the same at each point. A volume-average over the entire 
mesh approximates an ensemble average, carry out this volume average 
by a summation over the entire mesh and indicate the average of a by 
<a>. 


In Chapter I we spoke of the equivalence of Favre (mass weighted) 
and conventional Reynolds averaging in the homogeneous shear flow. This 
can be Illustrated rather simply as follows. 


Using the definitions in Chapter I, we write the total velocity as 
the sum of mean and fluctuating parts without specifying the type of 
average. Recall that the Favre average u^ is defined by u^^ = 
pu^/p and that the Favre and conventional fluctuating quantities, u^ 
and u|, are the differences between the unaveraged velocity and the 
respective mean. 


u. 






(4.3.1) 


^i = '^i + ("i " 


Using the definitions of the averages. 


PUf = pu^, pu^ = P Uj, 


pu^ = p u^ + pu|^ 


u^ = 0, pu^ = 0, 


= p Uj^ + p’u^ 


or 


u. 


= u, + 


P'u^ 




(4.3.2) 


(4.3.3) 
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The last term in (4.3.3) is the integral of the fluctuating momen- 
tum. Momentum is analytically conserved in a periodic domain, so that 
this term is a constant in time. The two velocity fields are then rela- 
ted by an additive constant. For convenience we choose the constant to 
be zero , which makes the Favre £ind conventional averages identical . 

4 .4 No nlinear Least-Squares Iteita-Flttlng 

In many places in this work we use a least-squares fit to data as a 
tool to study its behavior. Wi have 21 flow fields from which we take 
measurements. As discussed above, these flow fields span a range of 
shear, Mach, and Reynolds mmbers . this allows us to least-squares fit 
the 21 realizations of a particular quantity with a function that will 
show us the statistical variation of the quantity with these numbers. 

The fitting function that we use varies somewhat, depending on the 
quantity that we fit. In some cases we know what the behavior should be 
for the limits where one or more of the numbers becomes zero. The fit- 
ting functions reflect this known behavior. 

As an example we discuss the function 

f = d(^) (l+bM^)(Re^)^ (4.4.1) 

which has been used for a number of results . would use this func- 
tion for a quantity that would disappear when the shear rate became zero 
or infinity, depending on the sign of a — hence the appearance of 

/ CT 

( — ] . The Mach number term represents the first two terms in an expan- 
sion. This is a relatively standard analysis that may be found, for 
example, in Van Dyke (1975). It would be used for a quantity that would 
exist in an incompressible flow, but could be altered in compressible 
flow . 

The Reynolds number dependence is represented by the last term. In 
some cases we have arguments as to the behavior in the limits of zero 
and infinite Reynolds number. In this function, we have used a power 
law that disappears at zero or infinite Reynolds number, depending on 
the sign of c . 

The coefficient , d, is simply a scaling parameter. 
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As each function is used, we shall give a short explanation of why 
it has its form. The errors that are quoted are the r.m.s. deviations 
between the actual values and those predicted from the fitting function. 
These are then normalized on the actual values . They are percentage 
r.m.s. errors. 

4.5 Measures of Reynolds Stress Behavior 

life first discuss structural measures of the Reynolds stress tensor, 
< Pu^Uj >, in general and then seek changes that occur in the 
compressible flow. 

4.5.1 Shear-Stress Correlations 

The shear-stress correlation coefficient is defined 

2 2 1/2 

C = < uv > /[< u > < V > y ' (4.5.1) 

It is a measure of the strength of the turbulent shear stress and is 
experimentally known to be about C = 0 .5 in shear flows . Figures 4 .4a 
through 4.4h show the time development of this correlation coefficient. 
It will be noticed that the magnitude of the coefficient becmes quite 
large and then decreases throughout , the "good time” Indicating that the 
flows are still evolving throughout this period. There is too much 
scatter in the experimental data to say whether this trend is also 
observed in the laboratory. It is yet unclear in both numerical and 
laboratory experiments whether there is an asymptotic structural state 
to which these flows evolve. Evidence indicates that this may be so, 
but we are unable to carry the simulations far enough in time to deter- 
mine this. Further evidence will be given in support of the hypothesis 
that an as 3 miptotic state exists . 

The value of the coefficient during the good time is larger than 
experiment. Ihe average value over the 21 flow fields is C = 0.67, 
the standard deviation is 0.03. The corresponding value for the KJC 
(1977) flow is C = 0.47. The larger value that we calculate agrees 
with the simulations of Shlrani (1981) and of Rogallo (1979). 
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Suspecting that this difference may be a low Reynolds number 
effect, we least-squares fit the shear-stress coefficients with the 
fitting function. 

f = (1 + b(M)^)(Re^‘' (4.5.2) 

used this form because we expect that c will be zero if the shear 
number becomes zero, 1 .e . , there will be no shear stress if there is no 
mean velocity gradient. Wfe do not expect c to disappear at zero Mach 
number , and we look for Rejmolds number dependence as a power of Rey- 
nolds number. 

The values of the calculated constants in (4.5.2) indicate which of 
the nondlmensional parameters is most important in the variation of the 
shear stress coefficient. Kb find that the three constants a, b, c 
are essentially zero, indicating no dependence on the Reynolds number or 
the other parameters . Kb now suspect that these larger than experimen- 
tal values are due to the fact that the flow is still evolving. 


4.5.2 Principal Axis Measures 

The orientation of the coordinate system used in these simulations , 
although arbitrary, is the conventional choice. There is a coordinate 
system in which the stress tensor becomes diagonal. This is the prin- 
cipal axis coordinate system. The angle a between the two systems is 
defined as 


a 



-2 < puv > 

< pu^ > - < pv^ > 


(4.5.3) 


'rtiis measure of the stress tensor structure has been used by Corrsln. 
fe least-squares fit a by the fitting function, Eq . (4.5.2), in order 
to discover which parameter is most important in determining this rota- 
tion. 


Kb have again used (4.5.2) as a fitting function, but for a 
slightly different reason. At zero shear number a flow that has come 
to equilibrium will be isotropic. In this case (4.5.3) is indetermi- 
nate, because any orientation of the coordinate system is equivalent. 
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However, at infinite shear number the coordinate systems coincide. We 
expect a negative exponent on the shear number. We do not expect a to 
disappear with the Mach number, and we again look for Reynolds number 
effects with an exponent on Re;^. 

In Fig. 4.5, we plot the angle a against the estimate from the 
fitting function. The constants that were determined are shown in the 
upper left corner of the plot. If Eq . (4.5.2) were a perfect fit, the 
data would lie on a single 45° line radiating from the origin. Wfe see 
from the constants that the greatest dependence is on the shear number, 
as indicated by A = -0.40. For comparison we also plot the experimen- 
tal point of HOC (1977) on Fig. 4.l5 and note that it falls at one end 
of our correlation region, in good agreement with our results. 

In our simulations the principal axes of the stress tensor lie 
between 11 and 22° from the x-axis . The principal axis of the mean 
strain-rate tensor lies at 45°. Since eddy viscosity models force these 
two sets of axes to coincide, an eddy-viscosity model would not be 
appropriate in this flow. This conclusion also applies to one and two 
equation models, as defined by Reynolds (1976). These models calculate 
an eddy viscosity from the kinetic energy and a length scale. They 
cannot represent all the components of the Reynolds stress tensor at the 
same time and should not be used in flows in which more than one 
component of the stress tensor is Important. 


4.5.3 Principal Stress Ratio 

Wfe may also examine the ratio of the principal stresses . These 
stresses are related to the stresses in the unrotated corodinates by 

_ < pu^ > + < pv^ > 

” a,b 2 


/ < PU^ > - < p/ ^ ,2 


(4.5.4) 

2 

The transverse stress < pw > is also a principal stress and is not 
changed by the rotation. 

We calculate the ratio of the principal stresses in the X - Y, or 
shear plane, and least-squares fit these values with the function 
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f = (1 + bM^')(Re^)^ (4.5.5) 

life use (4.5.5) because we know the principal stress ratio is one in 
an equilibrium flow with a zero shear number. kb expect a positive 
exponent on (SL/q) . The arguments for the Mach and Reynolds numbers 
remain as before . 

In Fig. 4.6, we plot this ratio vs. an estimate from the fitting 
function (4.5.5) and show the values of the calculated constants in the 
upper left corner. Wfe again note a large dependence on the shear number 
indicated by a = 0.74. For comparison we also plot the HGC (1977) 
value . 

For a shear flow, Eqs. (4.5.2) and (4.5.5) could be used to check 
the effectiveness of a model in reproducing these measures of turbulent 
structure. 


4.5.4 Invariants of the Reynolds Stress Anisotropy Tensor 

As described by Lumley (1970, 1977, 1978, 1980), the invariants of 
the Reynolds stress anisotropy tensor can be used to characterize the 
Reynolds stresses. Ife first define the anisotropy tensor as 


ij 

ij 


\k ^ 


(4.5.6) 


This tensor is symmetric, traceless, has bounded maximum and minimum 
values, and vanishes when the stress tensor is isotropic. 


In Figs. 4.7a through 4.7h, we show the time development of the 
four nonzero components of this tensor in the shear flow. In all cases 
the diagonal components seem to approach an asymptotic state faster than 
the shear stress component, life cannot show that the shear flow comes to 
structural equilibrium (constant values of , because we are unable 
to carry the simulations further in time. However, the results strongly 
suggest that this is the case. In a later section we shall use this idea 
to derive a class of models for the pressure strain terms in the Rey- 
nolds stress equations . 
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The values of like those of ^ij* depend upon the coordi- 
nate system. Ive seek a measure of the stress tensor that is "Invariant" 
with the coordinates. Wh do this by contraction of various powers of 

and define 

I = bfi ; II = b^j b^j ; III = b^j^. b,^.j b^^^ ; etc. (4.5.7) 

From the Cayley- Hamilton theorem, as discussed by Lumley (1970), we 
find that I, II, and III are the only independent invariants of j • 
In addition, b^j is traceless, so 1=0. fe are left with two inde- 
pendent scalars to characterize the stress tensor. 

Figure 4.8 is an adaptation of Fig. 1 in Lumley and Newman (1977). 
In it we plot the values of II vs. Ill for the 21 simulated flow fields. 
Ife also plot the value for the HGC flow and show the limits within which 
all turbulence must lie . 

In seeking a key to the variation of I and II, we least-squares fit 
both of the invariants with the function (4.5.2) for the same reasons 
that we used in Section 4.5.1. Ihe results of this fit are shown in 
Table 4.3. Lfe again find the shear number to be the most important 
nondlmensional parameter in determining this measure of stress-tensor 
structure . 

We have consistently found that the shear number is the most impor- 
tant nondlmensional parameter in determining the structure of the stress 
tensor. These fits can be used as correlations in methods of predicting 
shear flows. In many flows the shear number is fixed, and it may be 
possible to base a model on this "structural similarity." 

4 .6 Direct Measures of Compressible Behavior 

As discussed in Chapter I, some of the original approaches to 
studying compressible turbulence involved decomposing the flow fields 
into parts and studying their interactions. Moyal (1951) decomposed the 
velocity spectra into "eddy turbulence" and "random noise" parts. For 
very low turbulence levels, he analytically calculated the Interaction 
of these parts to second order. 
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Numerical simulations allow us to perform these decompositions 
directly on the simulated flow fields and calculate their interactions 
without a limit on the strength of the turbulence. Wb carry Mayal's 
ideas further by using his spectral decomposition to define two real- 
space velocity fields. Wfe shall label the incompressible part (sole- 
noidal) u? and the compressible part (dilatation) uf* 

In Fourier space, a solenoldal (dilatation-free) velocity vector is 
perpendicular to the wave number vector. Ivb can decompose the Fourier 
transform of the velocity vector into the sum of a vector parallel with 
the wave number vector, u® and one perpendicular to it u|. Upon 
Inverse-Fourier transforming, we have the velocity fields u^ and 
u^. Operationally, we perform this decomposition as follows. fe define 

S,D_ S,9, / ) c. ^ \ 
u^ = gif (4.6.1) 

where is identified with the gradient of an undetermined scalar 
potential. Taking the divergence of (4.6.1), we find a Poisson equation 
for (j) . 


v2(j) = 




(4.6.2) 


Ihe solenoldal velocity field is calculated as the difference, u| = u^ 



Wfe can then decompose the Reynolds stress tensor into compressible 
(divergence) and incompressible (solenoldal) parts. Substituting the 
first equality in (4.6.1) into the definition of the Reynolds stress and 
dividing the density into mean and fluctuating parts , we have 


< PU^Uj > 



. S D ^ D S ^ 

< p(u^Uj + U^Uj + 


D D. . . . , . 

UiUj) > + < P U^Uj > 


R 


ij 


R 


D 

ij 


(4.6.3) 


where R^^ is the incompressible part and R^j is the compressible 
part of the stress . 


In Figs. 4.9a through 4.10h, we show the time development of the 
decomposed stress tensor for the 21 simulated flow fields. Only the 
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i = 2, j = 2 component of the stress tensor h as a significant com- 
pressible component . We believe that R 22 is related to acoustic 
noise, which is known to propagate preferentially in the Y or gradient 
direction away from a mixing layer. 

Wfe see similar behavior in the calculation of 3-D spectra of the 
decomposed components of the stress tensor . fe divide j by the 
density. In Fourier space the stress tensor becomes 

< u. u. > = < u.u. > + < (u.u. + u. u. + u, u.) > (4.6.5) 
Ij ij ''ijijij^ 

Incompr . Compressible 


Wfe have written this tensor as the sum of compressible and incompress- 
ible parts . 

The spectra of the terms in (4.6.5) are then calculated by the 
method described in Chapter III. The incompressible spectra contain 
information only from the solenoidal velocity field. The compressible 

q 

spectra contain the remainder. Wh designate these spectra as E£j 
and respectively. 

In Figs. 4.11a through 4.12h, we show the 3-D spectra of the dia- 
S D 

gonal components of and E^j at nondimenslonal time St = 6. The 

solid lines represent the spectra of the trace of spectrum of 

the trace of E^j is the standard 3D energy spectra, designated E(k) . 
It represents the distribution of the kinetic energy per unit mass over 
the range of turbulent length scales in the flow. Most of the informa- 
tion, as we have already seen, resides in the solenoidal component of 
the velocity field. This can also be seen by comparing the values on 

q 7^ 

the ordinates of the two sets of spectra, E^j and ^i j • 

The compressible energy spectra are dominated by the 1 = 2, j = 
2 components , especially at high wave number . Searching for a reason 
for this behavior, we formed the dynamic equations for pro- 

duction term appears to be responsible for the dominance of this compo- 
nent. So we write only this term 



... H — Skj^ kj^ u^,(u2 + 

(no summation) 
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(4.6.6) 



Note that this equation allows for production only when i = 2. 
IMs may be the mathematical reason for the dominance of the i = 2, j = 
2 term and it is probably associated with noise production. 

For completeness, we present the 3-D spectra of the shear stress as 
decomposed by (4.6.5). Again, most of the information is carried by the 
solenoidal part of the velocity field. These results are presented in 
Figs. 4.12a through 4.13h. 

Vfe have used a velocity decomposition based on Msyal's ideas. By 
dividing the velocity field into compressible and incompressible parts , 
we are able to decompose the Reynolds stress tensor . Our general con- 
clusion is that the stress tensor in compressible flow is very similar 
to its counterpart in incompressible flow. The only exception is the 
1 = 2, j = 2 component which is associated with noise production. 
These results are in keeping with Msrkovin's hypothesis that compressi- 
bility has little effect on the structure of the stress tensor. 

However, turbulence models still do not perform very well in high 
Mach number flows . suspect that the problem may be found in the 
modeling of the terms that contribute to the stress tensor. Perhaps 
these models cannot properly represent some compressible effect that 
occurs at higher Mach number, and, therefore, they do not produce the 
incompresslble-llke stress tensor that we have just examined. 

In the next section we discuss the dynamic equations for the Rey- 
nolds stresses in hopes of finding this problem. 

4 . 7 Reynolds Stress Equations 

The Reynolds stress equations are the time-dependent governing 
equations for the Reynolds stresses. As stated earlier, modeled ver- 
sions of these equations are solved in conjunction with the Navier- 
Stokes equations in the approach called stress-equation modeling. This 
more complicated type of modeling is introduced in hopes of capturing 
more of the dynamics of the Reynolds stresses and has found a great deal 
of favor in the literature. Although not yet an engineering tool, it 
shows promise in replacing lower-order models. 
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The stress equations are formed by taking moments of the Navler- 
Stokes equations. Due to the nonlinearity of the parent equations, the 
stress equations also contain correlations between fluctuating quan- 
tities. Unfortunately, these new correlations are of higher order in 
the velocity than the Reynolds stresses . If we are to use the Reynolds 
stress equations, we must model some of the terms that they contain. 
These modeled terms are generally not measurable by experiment — hence it 
is attractive to measure them in a simulation. 

The Reynolds stress equations for a homogeneous flow are written in 
Eq. (4.7.1). These equations incorporate the linear coordinate trans- 
formation described in Chapter II. The spatial derivatives that appear 
in (4.7.1) include the mesh metric St and are defined in 2.3.2. 

Ij- < pUjUj > = - < puju,^ > “j.k. ■ < > “l.k 






Pressure Strain 

(4.7.1) 

Homogeneous Dissipation 


TT P < u, i,(u. .+u . j) > Dilatation Dissipatl 
j K ,K. 1 , J 3,1 


on 


It • 

3t ij 


= ^3^ 


4) - (D^ + ) 

•^ij '‘"'ij "'ij'' 


(4.7.2) 


Labels for each term are included . 

As can be seen, the terms divide into three types: production, 
pressure strain, and dissipation. The convective terms do not appear in 
homogeneous flows , so we can concentrate on the homogeneous terms . 

Inspection of (4.7.1) shows that the production term, j » 
product of the Reynolds stress and the mean velocity gradient. These 
terms are responsible for amplifying the stresses and generally drive 
the flow away from isotropy. Since the stresses appear explicitly, 
along with the mean velocity gradient, these terms require no modeling. 
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The pressure strain term, involves correlations between the 
fluctuating pressure and the fluctuating strain rate tensor. Most at- 
tention in stress equation modeling is focused on these terms . They are 
responsible for redistributing the information among the components of 
the stress tensor. In the homogeneous shear flow, the pressure strain 
terms drain energy out of the i = 1, j = 1 component (where the only 
production occurs) and distributes it to the i = 2, j = 2 and i = 3, 
j = 3 terms, which have no direct sources of production. They also 
have a large effect on the 1 == 1, j = 2 or shear stress term. It is 
thought that at high Reynolds number this term is responsible for most 
of this redistribution. In compressible flows, this tensor is not 
traceless and therefore can have a net effect on the turbulent kinetic 
energy. The pressure strain terms must be modeled in a simulation. 

The homogeneous dissipation tensor, is familiar to incom- 
pressible turbulence modelers . At high Reynolds numbers it is thought 
to be nearly isotropic. However, it can have a redistributive character 
much like the pressure strain terms, if it is not. It is responsible 
for destroying the Reynolds stresses, and its trace dissipates turbulent 
kinetic energy. like the pressure strain terms, it must be modeled. 

The last term in (4.7.1) is the dilatation dissipation, j • This 
term occurs only in compressible flows. It is composed of the fluctu- 
ating dilatation and the fluctuating strain rate, and also requires 
modeling . We shall show that this term is small in comparison to the 
homogeneous dissipation term and that it may be neglected . 

Wfe shall discuss the behavior of the terms in the stress equations 
as calculated from the simulations, and then discuss some models used to 
approximate these terms . 

4.7.1 Time Behavior of the Reynolds Stresses 

We now discuss the time history of the Reynolds stress equations, 
as simulated in the compressible homogeneous shear flow. 

The contributing terms in the Reynolds stress equations (4.7.1) are 
plotted vs. the nondimensional time, St, in Figs. 4.14a through 4.17h. 
Ml terms have been normalized by the shear rate times the trace of the 
stress tensor. 
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In this geometry, the production terms operate only in the i = 1, 
j = 1 and i = 1, j = 2 equations, i.e., there is production only in the 
the equations for < pu > and < puv >. This is seen by examining the 
figures . These terms increase the magnitude of < pu'* > and < puv > . 
During the good time of each simulation, each of these terms has become 
relatively constant . 

The pressure strain terms show the expected behavior in the shear 

2 

flow. The term in the < pu > equation acts to drain energy from 

2 

< pu > . The pressure strain terms are the only terms that contribute 
2 2 

positively to pv and pw . In the shear stress equation, the pressure 
strain term acts to destroy the stress. 

The homogeneous dissipation terms act to destroy the stresses in 
all four equations. Note that the dilatation dissipation is virtually 
zero in all cases. In comparison to the homogeneous dissipation, it is 
very small and, therefore, we neglect it in the rest of the analysis. 

In Figs. 4.18a through 4.18h, we show the time development of terms 
in the turbulent kinetic energy equation. This equation results from 
taking the sum of the diagonal Reynolds stress equations . fe show these 
data to point out the dissipative behavior of the pressure dilatation 
term. In an incompressible flow, the pressure strain tensor is trace- 
less because the velocity field is dilatation-free. However, in com- 
pressible flow, the pressure dilatation term may have a net effect. Our 
simulations show this term to be dissipative of turbulent kinetic 
energy. This effect is not represented in turbulence models derived 
from the incompressible equations, but we feel that it should be. 

4.7.2 Dissipation Anisotropy 

In high Reynolds number flows, the homogeneous dissipation tensor 
D^j (hereafter referred to as just the dissipation tensor, D^^) is 
thought to be isotropic, tfe investigate this by forming the dissipation 
anisotropy tensor in the same way that we formed the stress anisotropy 
tensor . 

*^ij ° (4.7.3) 
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If the dissipation is isotropic, then would be zero. 

In Fig. 4.19, we plot d^^j vs. the stress anisotropy j * ^ 

find that they are nearly proportional. '[his is a very surprising 
result, which we first attributed to the low Reynolds number of our 
simulations. It does not seem to be spurious in that it is also seen in 
the simulations of Shirani (1981) and of Rogallo (1979) . have a 

Taylor microscale Reynolds number range of about 18-120, so we thought 
that we could search for Reynolds number dependence by studying the 
Invariants of d^|^ j . These invariants should disappear as the Reynolds 
number increases . 

Ife form the invariants in the same was as for the stress tensor and 
least-square fit these scalars by the fitting function. 

f = d (1 + bM^)(Re^)^ (4.7.4) 

Ife chose this function to look for power-law dependence on both shear 
and Reynolds numbers and because we do not expect the invariants to 
disappear with Mach number. 

Ivfe expected to find a positive exponent, a, on the shear number, 
and a negative exponent, c, on the Reynolds number. Ife were surprised 
to find a positive exponent on the Reynolds number. The results are 
shown in Table 4.4. This indicates that there is no tendency for the 
dissipation to become isotropic at high Reynolds number. This is a very 
Important result that we believe may be characteristic of shear flows . 
Ihis fact does not seem to be generally recognized and is very important 
in turbulence modeling. It plays an important role in some of the 

models that we subsequently evaluate. 

4 .8 Character of the Pressure 

Before discussing the pressure strain terms, we must look at tVie 
character of the pressure field Itself. 

In compressible flows the pressure field obeys a hyperbolic equa- 
tion. It has a wavelike character and information travels at a finite 
speed. Incompressibility results from assuming that this speed of 
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propagation is very large in comparison to convective velocities. This 
is reflected mathematically by the elliptic nature of the pressure in 
incompressible flows, in which the pressure field adjusts instantane- 
ously everywhere to changes in the velocity field. 

This can be shown by taking the divergence of the incompressible 
momentum equations to form a Poisson equation for the pressure. 


1 V^p 

P 


u. .u . , - u. .U . , 
i.J J,i i.J J,i 


(4.8.1) 


The Poisson equation is elliptic. In incompressible flows, it is solved 
in conjunction with the momentum equations to enforce the condition of 
zero dilatation or mass conservation. Equation (4.8.1) shows the pres- 
sure to be completely determined by the velocity field. It is not 
really an Independent variable . In compressible flows the pressure is 
determined from the thermal energy equation. It is truly a flow vari- 
able and represents an additional degree of freedom. 

\k may show the connection between the two ways of calculating the 
pressure by taking the divergence of the momentum equation (2.3.4). 


2 

V^P + 


3t 


(Pu^)^i 


= - (pu.u.) ,, - (pu, ) ,U, 




i ,3 J,i 


(4.8.2) 


Using the conservation of mass equation (2.3.3), (4.8.2) becomes 

2 

V^P - — ^ p = g, + g, (4.8.3) 

3t 

where gj, and g 2 represent the right-hand side terms of (4.8.2). 
Using the definition of the speed of sound, we may write (4.8.3) as 

V^P - = ®1 (4.8.4) 

3t^ c 

which is a nonlinear. Inhomogeneous wave equation. This illustrates the 
wavelike or hyperbolic nature of the pressure, in contrast to the ellip- 
tic behavior in incompressible flows. 

As the sound speed, c, becomes very large, the time derivative 
term in (4.8.4) becomes small, and in the limit we are left with the 
Poisson equation (4.8.1). 
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It is emphasized that all pressure fields in these simulations are 
calculated from the full compressible Navier-Stokes equations . 

Turbulence modelers have exploited the fact that the pressure field 
is the solution of a Poisson equation. In the next several sections we 
apply this reasoning to the compressible flow in order to relate the 
pressure strain terms to commonly used models that were derived from the 
incompressible equations . 


4 .9 Poisson Decomposition of the Pressure Field 

Existing pressure strain models have been derived from the incom- 
pressible flow equations . Historically, modelers have decomposed the 
pressure strain tensor into two parts, in order to construct models 
individually for each part. The Poisson equation (4.8.1) is used to 
divide the pressure field into two parts, which are then combined with 
the fluctuating strain rate to form two pressure strain tensors . 

In order to relate existing models to our simulated flow fields, we 
follow the same process in the compressible shear flow. Remembering 
that the mathematical nature of the flow field has changed, we decompose 
tVie actual pressure strain field in the manner that we shall describe. 
Ttiis can be done only in a numerical simulation. Ihis is the first time 
that these terms have been directly calculated. 

Wfe take our direction from (4.8.3). If we define 


g2 

go = —9 P (4.9.1) 

3t^ 


we then rewrite (4.8.3) as: 


V^P = 


?1 + 83 + 83 


= -(Pu^u.)^^. - (pup^.U.^^+ 


(4.9.2) 


Rotta 


Fast 


8t 

Compressible 


Th;ls is a Poisson equation for the pressure, with g 3 written as a 
source term representing the compressible contribution to the pressure 
field. Equation (4.9.2) does not reflect the hyperbolic behavior of the 
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pressure, but we use it here to make a connection with models commonly 
proposed for Incompressible flow. 

note that (4.9.2) is linear in the pressure. Therefore, the 
pressure field can be regarded as the sum of three fields derived from 
each of the three source terms on the right-hand side of (4.9.2). Ife 
write; 



(4.9.3) 


Combining each pressure with the fluctuating strain rate, we have 


k‘*’ij 


= < 


V"i,j-^ '^j,i> > 


(4.9.4) 


where we have again used the symbol < > to denote volume average, tfe 
now have three distinct pressure strain tensors that sum to the total 
tensor. Vfe identify each piece in keeping with its pressure term in 
(4.9.2) . 

is the "Rotta" term involving a pressure resulting from the 
turbulence interacting with itself, originally identified by Rotta in 
1951. 2*^1 j mean-fluctuating term involving a pressure resulting 

from interactions between the turbulence and the mean velocity. We 
hereafter call this term the "Fast" term (Lumley, 1978) . 
call the compressible component, as it accounts for the wavelike behav- 
ior and does not exist in the incompressible flow. Other decompositions 
are possible. The inclusion of the fluctuating density alters the defi- 
nition of the Rotta and Fast terms slightly; however, their essential 
character is preserved in making this analogy. 


4 .10 Magnitude of the Pressure-Strain Terms 

It is customary to propose models for each part of the pressure- 
strain tensor separately, and then use them in combination to replace 
the total pressure strain tensor . The performance of the model is then 
judged indirectly by its effect on measurable quantities, such as the 
stresses or the mean velocity. 
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Wb find ourselves in the fortunate position of being able to 
directly calculate these pressure-strain terms and compare models to 
them Immediately. decompose these terms in the manner suggested by 
Eqs. (4.9.2) through (4.9.4). 

Having expected the compressible part to be small in comparison to 
the Rotta and Fast parts, we are surprised to find it of the same order 
of magnitude as the others . As a measure of the magnitude of this part 
of the tensor, we form the second invariant and divide it by the sxim of 
the second invariants of the Rotta and Fast parts. Wfe performed a 
least-squares fit to this quantity, to look for the important parameters 
in its development. In Fig. 4.20, we show this ratio of invariants 
plotted against the least-squares fitting function. The fitting func- 
tion and the determined constants are also shown. There is a positive 
dependence on both the shear number and on the fluctuating Mach number. 
Recall from Section 4.1 that the product of these two parameters is the 
shear Mach number. It seems as though the magnitude of the compressible 
part of the pressure-strain term probably scales on this "external flow- 
like" Mach number. In these simulations, the magnitude of this term 
seems to vary almost linearly with the shear Mach number, as indicated 
by the exponents a and b. This term will become important in mixing 
layers and boundary layers in which the Mach number difference across a 
large eddy becomes greater than one . 

It will be useful to describe the procedure used to verify this 
result. The procedure used involves the direct calculation of the first 
three terms in (4.9.2). The fourth, the compressible source term, is 
calculated by difference. Equation (4.9.2) is analytically satisfied by 
the flow field . It is numerically satisfied if the same method of spa- 
tial derivatives is used for it as for the original simulation, kb have 
used pseudo-spectral derivatives in all cases as described earlier, for 
numerical consistency. 

The "compressible" part of the pressure will be very small in a 
low-Mach-number isotropic flow. kb know that, at low Mach numbers, this 
flow is virtually incompressible. In the simulation labeled IH64A, we 
calculated the ratio of the pressure fluctuations for a flow field at 
the end of the computation, when the fluctuating Mach number was less 
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than 0.03. The second portion of the pressure does not exist in this 
flow, as there is no mean velocity gradient, fe found the ratio of the 
squares of the fluctuating pressures to be 

= 0.0056 (4.10.1) 

This is quite small, consistent with our hypothesis. 

The third source term was calculated from the divergence of the 
momentum equations. As a check, we evaluated g^ in an independent 
way. From (4.8.2) and (4.8.3), it can.be seen that this source term is 
the second time derivative of the density. The code is run two time 
steps to produce three consecutive flow fields, and 9^/3t^p was then 
approximated by a three-point central difference formula. Ite then found 
this ratio of pressure fluctuation to be 

Pj^/P{^ = 0.0050 (4.10.2) 


iin good agreement with (4.10.1). The magnitude of (4.10.2) is slightly 
less than (4.10.2), as might be expected, because of the high wave 
number dissipation inherent in the three-point formula. 

This confirms that our hypothesis and method of calculating the 
pressure is correct. At very low Mach number, the compressible part of 
the pressure disappears and, with it, its contribution to the pressure- 
strain terms . 

In order to find out what physical mechanism was responsible for 
the size of P^ , we applied the chain rule to the source term g^ and 
calculated four pressure fields that comprise P 3 . 






This was done in the simulation HS64D at St = 6.0. again formed the 
fluctuating pressures and found that the last term in (4.10.3) was re- 
sponsible for more than 95% of the P^ pressure field. The other three 
terms are insignificant . This Indicates that the third or compressible 
part of the pressure field is a result of the velocity dilatation in the 
flow field and not a result of variations in density. 


66 



4.11 Ifodellng of the Decomposed Pressure-Strain Terms 

In this section we discuss and evaluate models for the three compo- 
nents of the decomposed pressure-strain tensor. 

4.11.1 Rotta Pressure-Strain Term 

The Rotta pressure-strain term, 1‘t’ij* first identified by 

Rotta (1951). The pressure field associated with this tensor results 
from Interactions of the turbulent velocity field with itself . It may 
therefore exist in a flow which has no mean velocity. In a flow with no 
turbulence production, this term is thought to be responsible for the 
experimentally observed return to isotropy of the Reynolds stresses . 

In order to model this part of the pressure strain tensor, Rotta 
related the to the stress anisotropy, by reasoning that 

the effect of the model should stop when isotropy is reached. He pro- 
posed the simplest model that is linear in 

I't’ij “ ^^ij (4.11.1) 

where e is the isotropic part of the dissipation defined as 

e = y (4.11.2) 

Although Rotta proposed this model solely as a replacement for we 

shall show that there are other terms that should be combined with the 
Rotta term and modeled this way. can also make an estimate of the 

value of C2^ • 

Consider the time- dependent equations for bj|^ j , the Reynolds 

stress anisotropy tensor, in a homogeneous flow with no turbulence pro- 
duction (no mean velocity gradient) . Ife do this to isolate the effects 
of the Rotta term only. In this relaxing flow, there is no fast 
pressure-strain term, 2^ij> pressure-strain tensor remains the 

sum of the Rotta part and the compressible part. laray (1980), however, 
has shown that the velocity dilatation quickly becomes insignificant in 
a relaxing, compressible, homogeneous flow and therefore the pressure- 
strain tensor, in this flow, is almost entirely due to the Rotta term. 
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In addition, the trace of is small, so its dissipation may be 

neglected for this case. 

Manipulation of the Reynolds stress equations produces the Reynolds 
stress anisotropy equation. 

I- b, . - -=^r, - 2Ed, , 1 + -^ Eb, . (4.11.3) 

3t Ij Uj 13 

where we have used d^j for the dissipation anisotropy, defined by 

(4.7.2) and assumed that is traceless. 

In (4.11.3) the term in square brackets can alter the relative mag- 
nitudes of the elements of l«e., it has redistributive behavior 

and is responsible for the return to isotropy. The last term cannot do 
this. This indicates that the entire bracketed term should be replaced 
by a model rather than just the pressure-strain term. 

If we use the Rotta model, (4.11.1), to replace fl'l’ij “ 2edj^j] , 
we find 

fcbjj = (2 -Ei) b^j (4.11.4) 

For b^j to relax to isotropy, Cj^ must be greater than two. 

Thus it makes sense to use the Rotta term in place of the term in 
square brackets. The constant will be evaluated for both applications. 

life are able to directly calculate the constant Cj^ from our simu- 
lated flow fields . All three pressure-strain tensors exist in the homo- 
geneous shear flow, so we calculate the Rotta pressure-strain term by 

(4.9.2) through (4.9.4). 

In Table 4.5 we show the value of the constant C]^ for each of the 
four nonzero Reynolds stress equations. Note that, for both applica- 
tions, Cj^ is not a constant. It varies among the Reynolds stress 
equations and has a spread of values for each equation alone. Ife fit 
the constants, Cj^ , calculated from each flow fields with the following 
two equations ; 

f = 2 + 1 (1 + bM^)(Re^)‘^ (4.11.5a) 
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f 


d 


(4.11.5b) 


(1 + bM^)(Re^)‘' 

in hopes of determining which nondlmensional parameter varies with the 
constant. The 2 appears in the first fitting function because cj^ 
should approach 2 (no return to Isotropy) at zero Reynolds number, at 
least for the case where the dissipation anisotropy is included. This 
argument is discussed in Lumley and Newman (1977) and motivated the 
first form. Ivfe were unable to find a suitably converged set of con- 
stants for this equation, which led us to the second form in (4.11.5), 
which we found to be a better fit to the data . These results are pre- 
sented in Table 4.5. 

From the values of the fitted constants in Table 4.5 it can be seen 
that the shear number is important in determining the variation of Cj^ . 
It seems as though the mean velocity gradient has an indirect effect on 
the structure of this term. 

The values of Table 4.5 are shown graphically in Figs. 4.21 and 
4.22. Figures 4.21 and 4.22 are actually four plots in one. We deter- 
mined the constants in Eq . (4.11.5b) by least squares. Wfe then used 
these constants in (4.11.5b) to make an estimate of c^^ in each flow 
field. Wfe divided the predicted value by the actual value and plot this 
ratio on the ordinate. If the least squares were perfect, this value 
would be one. For clarity, we shift the origin for each plot by 
successively adding two to the results from each of the four stress 
equations. The results for the i = 1, j = 1 equation appear at the 
top of the plot beside the number 7 . On the abscissa we plot the 
corresponding values of Cj^ . Wh may then visually compare the range in 
values of Cj^ and the accuracy of the fitting function for the four 
stress equations. 

If the model were a perfect representation of the terms it re- 
places, the figures would have only four points, each immediately aboye 
the others, at a common value of cj^ . However, the horizontal spread of 
values shows the variation of cj^ in each equation, and the vertical 
spread shows the adequacy of the least-squares fit . 

The values of Cj^ for the case where the dissipation anisotropy is 
included are, within the variance, greater than two, consistent with the 
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requirement for return to Isotropy. For the simple Rotta model, this Is 
not true. (Note that, when Is evaluated from experiment, the dis- 
sipation anisotropy Is usually Included with the Rotta term by default, 
because the dissipation Is assumed to be Isotropic. Thus the resulting 
models are fairly accurate , even though some of the assumptions are not 
correct .) 

In the relaxing, homogeneous simulations of Rogallo, Cj^ Is cal- 
culated to be constant at a value of about 2.37. The values of Cj^ 
calculated by Schumann and Patterson (1975) are of the same order of 
magnitude as we calculate for the simple model. They do not Include the 
dissipation anisotropy In their analysis. The Incompressible, homogene- 
ous shear simulations of Shlranl (1981) also show a variation of Cj^ 
when the Rotta term Is calculated directly. It seems as though the 
Rotta model breaks down In the shear flow. 

Bearing this In mind, we have looked at some arguments about the 
variation of c^. Various workers have recognized that Cj^ Is not a 
constant and have proposed arguments as to how It should vary. Lumley 
and Newman (1977) argue that c^ should be a function of the Reynolds 
number and the two nonzero Invariants of the stress anisotropy tensor. 
They propose the equation 

/T, 1/2 

-a/Re ^2 / i . 
c^ = 2 + e ^ “TTI ^ c(dIII - II)) (|-+ 3III + II) 

(4.11.6) 

They chose the form of this function by determining how Cj^ should 
behave In the limits of zero and Infinite Re and anisotropy 
Invariants. The equation Is simply an Interpolation between these 
limits. fcfe performed a least-squares fit to the four constants in 
(4.11.6) and present these data along with error estimates In Table 
4.6. life were unable to find a converged set of constants for the 1 = 
3, j = 3 equation, but found reasonable values for the other three. 
The values that Lumley proposes are also Included . The quantitative 
agreement with his constants Is not good, but they are all of the same 
order of magnitude . 

Chung and Alrian (1979) use a similar argument to propose the 
fitting equations : 
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(4.11.7) 




2 + 


PTRe^) 


(-bll) 

all e (1 - 


+ , 111 ) 


u - 1 + f 

ij 

Ufe; were unable to find a converged set of constants for this equation. 

A simpler suggestion by Lumley and Newman (1977) does not depend on 
the stress invariants . 


C]^ = 2 + A Re^^/^ (4.11.8) 

Table 4.7 shows the values of A that we calculate from the shear flow. 
Note the variation and especially the change of sign. 

have to conclude that the Rotta model (4.11.1) is inadequate to 
represent the Rotta and dissipation anisotropy terms, when analyzed by 
this method. The constant Cj^ is a function of i and j and perhaps 
should be a tensor. As previously stated, the Rotta model is linear in 
bj[^ j . Perhaps higher-order terms in b^ ^ should be investigated that 
would represent this variation with a scalar constant. Ife suspect that 
the Poisson decomposition is at fault and that perhaps the Rotta model 
is not as bad as this analysis shows. Wfe shall return to this model 
when we attempt to model the entire pressure-strain tensor . 

4.11.2 Fast Pressure-Strain Term 

The Fast portion of the pressure-strain tensor was recognized by 
Rotta (1951), but its name is due to Lumley (1978). It arises from 
interactions between the turbulence field and the mean velocity gradi- 
ent . Many workers have proposed models for this term which have grown 
in sophistication over the years . Madern models are generally construc- 
ted from combinations of there is some controversy as to 
whether these models should be linear or of higher order in bj|^.. 
(Lumley, 1978). 

The reasoning that leads to the use of bjj^j in these models was 
given by Lumley (1974). Assuming incompressible homogeneous flow, we 
may write the Fourier transform of the fast part of the pressure as 
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(4.11.9) 


/V -ik.u. 

p = u 

'2 ^2 “j,i 

The pressure-strain term is written as an integral in Fourier space . 

< P^Cu^^j+Uj^^) > = ly* jlkj(P* u^+Pu*) + ik^(P*Uj+Pu*)jdk (4.11.10) 


Inserting (4.11.9) into (4.11.10) and remembering that U- ^ is a con- 
stant, we find 

Ic • A Ic 




dk 


“ ^^k^ijkJl (4.11.11) 

Equation (4.11.11) is an exact representation of the fast term for an 
incompressible homogeneous flow. It shows that the fast term is com- 
posed of a fourth-rank tensor multiplied by the mean-velocity gradient . 

Models for this term approximate the fourth-rank tensor because we 
cannot directly evaluate the integral. Note that the integrand is com- 
posed of Fourier-space representations of the Reynolds stress . This 
justifies constructing models for the tensor from combinations of the 

"ij- 

The most sophisticated model in use is linear in b^j . It has been 
derived in several forms by various workers, but all are essentially the 
same. Ife shall discuss the form found in Reynolds (1976). The tensor 
^ljk£ approximated by linear combinations of each term 
multiplied by its own undetermined constant. Constraints implied by 
(4.11.11) eliminate all but one of these constants, which must then be 
determined empirically. Mien combined with the mean-velocity gradient. 


U 


il,k« 




the model is 


(f * IT Vkl] 


(4.11.12) 


4 /5 

5 
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' 2 <"l,J - 


Itiis is exactly the form used by Hanjalic and Launder (1974) . 

For the homogeneous shear flow, the four nonzero terms are: 

2‘^U = " (t5 ^ 

2:‘*’22 “ il \ l )^ 

(4.11.13) 

2‘**33 ° (“ T 

2*^12 " 5 “ (l5 ^1 j) ^ ^ “ (it ^ " i) ^ ^ ^'^2 ^ 

As with the Rotta term, we can directly calculate the value of the 
constant Aj^ from our simulations. Table 4.8 shows the value of Aj 
for the four nonzero Reynolds stress equations. It can be seen that the 
average value of Aj^ varies among the four equations , but that the 
variance of Aj^ is quite small relative to the average for three of the 
four. The values calculated for the R 23 term have some scatter due to 
the small value of 1 * 33 , and are less reliable, life present them for 
completeness . 

To look for variations of the constant with the nondimensional 
parameters, we have least-squares fit the constants Aj^ from each of 
the four equations with the function 

= d (1 + bM^)(Re^)‘^ (4.11.14) 

Ife have chosen this form rather arbitrarily to look for power-law 
variations of Aj^ with shear and Reynolds numbers and because we do not 
expect Aj^ to disappear with Mach number. There have been arguments 
proposed that this "constant" should vary with the stress anisotropy 
also (Luialey, 1978). There is no clear indication ' of how it should 
behave in the limits of zero and Infinity for the shear and Reynolds 
number and for the invariants , and therefore no function has been 
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proposed, life justify (4.11.14) solely by the accuracy of the fit shown 
by the small rms error in Table 4.8. Figure 4.23 presents the con- 
stant in the same manner used to present the Rotta model constant. 

The constants are quite constant for each equation, leading us to 
suspect that much of the Important physics is captured by this model . 
However, the average value still varies from equation to equation. 
Improvements still need to be made, perhaps along the lines suggested by 
Lumley (1978), where he discussed Schumann's (1977) concept of realiz- 
ability . 

A simplified version of this model has found favor among some 
workers. It is known as the Gibson-Launder model, and is 


where 



< pu^uj^ > 'Jj.k - < > ^i,k 


(4.11.15) 


and 





li 


where ^ is the rate of turbulence energy production per unit volume. 
For the homogeneous shear flow, 

2*11 ■ ~ 

2*22 ■ 1*33 -1*2^ (4.11.16) 

2*^12 ~ ^ ^ ^“2 ^ ® 


The values of A 2 are presented in Table 4 .9 and graphically in 
Fig. 4.24. This model does not perform as well as the more general 
fourth-rank tensor model from which it is derived. One problem is the 
equality of the model for 2"^22 2'i’33* Ihere is no reason why 

these terms should be equal, and there is evidence that they are not 
(Champagne et al . (1970) and Harris et al. (1977), but the model is 
incapable of representing this. It seems as though some essential 
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physics has been left out by the simplification, and we would therefore 
prefer the more general model (4.11.12) over the Gibson-Launder model 
(4.11.15) . 

In hopes of Improving the fourth-rank tensor model, we have added 
terms involving the dissipation anisotropy. This addition is 

motivated by our experience with the Rotta term and also by the appear- 
ance of the integrand in (4.11.11). The wave numbers multiplying the 

velocities represent spatial derivatives and appear similar to the 
Fourier space representation of the dissipation tensor. 

Ife added terms linear in the dissipation anisotropy d^j into the 
model (4.11.12). Wb applyed the kinematic constraints discussed in Rey- 
nolds (1976) and also in Lumley (1978). Writing the model for the shear 
flow , we have : 

2*11 ■ (it ^ + t) <<^2 

2^22 " ”(t \ 

A 2 . /p^ °12 2 

2^3 = 5 C 2 spq 

2 ‘*’12 “ (1 + Aj^) spq2 - yT T ^ ^ ^'^l ^ 

- XY A^ - I s < PU 2 > + I C 2 + 2 C 3 spq^ (4.11.17) 

/ 0 ^ \ 2 ^ ”^22 2 

where is the dissipation tensor of EQ • (4.7.2)). Note the 

appearance of three undetermined constants , A 2 ^ , C 2 , C 3 . To determine 

these constants we did a least-squares fit to 2‘*‘ij* using (4.11.17). 
Title results of this fit appear in Thble 4.10. 

The constant A 2 ^ which multiplies the same terms as in the origi- 
nal model (4.11.12) varies much less from equation to equation than it 
does in the case where d^j is not included. This is an improvement. 
Hawever, C 2 and C 3 do not appear to be constants at all (c 2 and 
C 3 were determined simultaneously in the Rj^]^ and R 22 equations and 
are therefore equal) . Wh therefore do not expect this model to be an 
improvemtint over (4.11.12). 
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There are serious objections to the form that the fourth-rank ten- 
sor model takes (Leslie (1980) and Lumley (1978)), but In the absence of 
any demonstrably better model we would use (4.11,12) with carefully cho- 
sen Individual constants for each Reynolds stress equation. 

4.11.3 Compressible Pressure-Strain Terms 

The pressure strain terms associated with the compressible part of 
the pressure can be thought of as a deviation of the compressible terms 
from the equivalent Incompressible terms . 

In seeking a clue as to how they behave, we plotted the ratios of 
the Individual terms In to the corresponding terms In the Rotta 
and fast parts of the pressure-strain tensor. Ivfe noted that 3‘l^ij is 
similar to the fast tensor 2‘^'lj* ratio of these two terms seems to 
be a function of the shear Mach number SL/c, as shown In Figs. 4.25a 
through 4.25d, where we have plotted the ratio 3<l'ij/2‘)'i j four 
nonzero components vs. SL/c. This trend may be seen In all except the 
1=2, j = 2 term, where we see too much scatter to draw conclusions. 
Vt are led to suspect that the models for the fast term may be appli- 
cable to the compressible term. 

Vfe applied the most successful model (the general fourth-rank ten- 
sor model) for the fast pressure-strain terms to the compressible term. 
The results are shown In Table 4.11 and Fig. 4.26. By comparison with 
Table 4.4 and Fig. 4.23, It can be seen that the compressible terms 
behave In much the same way as the fast terms . 

Since these two tensors are quite similar, we thought that their 
sum should be modeled. implying the fourth-rank tensor model to 2*^1 j 
+ o<j)jj, we find the results shown In Table 4.12 and Fig. 4.27. The 
average value of does not vary as much among the equations as when 
applied to the compressible terms alone. 

Vfe evaluated several other models for these terms , Including the 
Rotta model and the fourth-rank tensor model with the dissipation anlso- 
tropy terms, but found none that performed as well as the original 
fourth-rank tensor model. 
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4.11.4 Conclusions for Madeling of the Decomposed Pressure Strain 
Tensor 

Using the simulated flow fields, we decomposed the pressure field 
by means of a Poisson equation. This is a standard decomposition used 
in Incompressible flow modeling and is used to give insight into con- 
structing models for the pressure-strain terns . Vfe carried out this 
decomposition in the compressible flow to relate existing pressure- 
strain models to the "exact" terms. 

Ife found three components instead of the two found in incompress- 
ible flows . The third component we labeled the compressible pressure- 
strain terms and regarded it as a deviation from incompressible flow. 

Vfe evaluated models against the exact terms, but found variations 
in the values of the constants. On the basis of these calculations, we 
would recommend the use of the Rotta model to replace both the Rotta 
term and the dissipation anisotropy. Vie recommend the use of the 
fourth-rank tensor model to replace the sum of the fast and the com- 
pressible pressure-strain terms . 

However, we have reservations about the validity of this Poisson 
decomposition in a compressible flow. As discussed in Section 4.8, the 
pressure is a hyperbolic quantity, not elliptic as in the incompressible 
flow. A Poisson equation introduces an elliptic character into the 
pressure and is perhaps Incorrect in a compressible flow. This leads us 
to two approaches described in the next two sections, where we first 
look into decomposition of the pressure by means of a wave operator, kb 
then evaluate several models for the entire pressure-strain tensor with- 
out decomposition. 


4 .12 khve-Operator Decomposition 


Because the Poisson decomposition is mathematically incorrect in a 
compressible flow, we have Investigated the use of a wave operator. We 
rewrite Eq . (4.8.3) 


V^P 


9^ P 
9t^ 


= 81+89 


(4.12.1) 
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where and §2 are the Rotta and fast source terms defined in 

Section 4.8. Ideally, we would like to calculate two pressure fields 
corresponding to the source terms gj^ and g 2 • Ife would then combine 
them with the fluctuating strain-rate tensor in the manner of (4.11.10), 
to form a Fourier space representation of two pressure-strain tensors. 
Equation (4.12.1) is a nonlinear wave equation which we cannot directly 
solve. However, we can retain the hyperbolic character of (4.12.1) by 
approximating the sound speed c as a constant, leaving a linear wave 
equation. 

Setting c = c^, a constant, and Fourier transforming (4.12.1), we 

find : 

^2 A, fy yy A A A 

P + (k c^) P = g, + g, (4.12.2) 

gt2 o 12 

For each source term g^^, the general solution of (4.10.2) is 

P^(k,t) = P^e ° + ?^e ° sin(t-t’)g^(t’) dt' (4.12.3) 

o*' 

Vfe are unable to explicitly evaluate the integral in (4.12.3). It rep- 
resents history effects that would appear in the pressure-strain terms . 
This is to be expected in that the pressure is truly a flow variable 
that develops according to its own dynamic equation. It is not a result 
of a kinematic constraint , as it is in incompresible flows and therefore 
is not completely determined by the instantaneous velocity field. 

It seems that models developed from this decomposition must have 
history effects built into them. This appears prohibitively expensive 
for useful models in that the entire history of the simulation must be 
preserved to construct the pressure-strain model. It is possible that 
models could be constructed on this basis, but we abandoned it for lack 
of time. 

We were therefore led to look elsewhere . In the next section we 
discuss attempts to model the entire pressure-strain tensor without any 
decomposition. 
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4.13 Modeling of the Entire Pressure-Strain Tensor 

Vfe have discussed our reservations about the use of the Poisson de- 
composition for the pressure-strain tensor in a compressible flow. Ub 
then described a decomposition based on a wave operator that appears 
impractical to use. 

Wa are led to models for the entire tensor without a decomposition. 
The models we evaluated follow two different lines of thought. Combina- 
tions of existing models were evaluated and their constants were simul- 
taneously determined. Vfe also constructed and evaluated models based on 
a structural similarity assumption. First we discuss combinations of 
existing models. Ivfe then discuss models based on a structural similar- 
ity concept. 

4.13.1 Combinations of Existing Models 

In a previous section we evaluated models against the exact parts 
of the decomposed pressure-strain tensor. Because of our reservations 
about the decomposition, we reevaluated the constants in two sets of 
models by sumultaneously determining them. Following our tentative 
recommendations at the end of Section 4.11, we model the sum of the 
total pressure-strain tensor and the dissipation anisotropy by the sum 
of the Rotta model and a fast model. 

/■A - ^ Fast model , Rotta model 

'■**’ij ij'' with constant with constant ^ ‘ ' 

The two constants are then determined by a least-squares fit using the 
models as a fitting function. W; did this hoping that any artificiality 
associated with the decomposition could be eliminated. 

Both the Gibson-Launder model and the fourth-rank tensor model were 
evaluated. The results are shown in Table 4.13 where we have written 
the values of the two constants and the normalized RMS error associated 
with the fit. 

The i = 3, j = 3 term has a large amount of uncertainty associa- 
ted with it because of the small and erratic value of the stress aniso- 
tropy in this term. kfe notice, however, that the three remaining terms 
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are better approximated by the fourth-rank tensor model. (There is less 
variation of the constants from equation to equation.) Comparison with 
Tables 4.4 and 4.5, where we have written the individually determined 
constants, shows less variation in their values fe take this to confirm 
our suspicion that a Poisson decomposition should not be used in the 
compressible flow. 

The combination of fourth-rank tensor and Rotta models is still the 
best one we have evaluated. Our recommendation of Section 4.11 still 
holds . 


4.13.2 Structural Similarity Models 

Wien analyzing their time history, we noticed that the elements of 
the stress anisotropy tensor seemed to approach as 3 miptotic values. This 
was particularly true for the diagonal terms, but less so for the shear 
stress. It was thought that we could exploit this fact to create a 
pressure strain model. 


As discussed earlier, it is by no means certain that the shear flow 

g 

comes to structural equilibrium b^^ = 0). However, the time rate of 

change of b^j becomes small during the "good time" of the simulations, 
and we use this to justify the zeroing of (9/9t)bj^j. This is an as- 
sumption similar to that used in Rodi's (1976) algebraic stress model. 


Wfe derive the time-dependent equations for b^j from the Reynolds 
stress equations (4.7.1) and (4.7.2). ^plying the chain rule to the 
time derivative of the definition of 


i_b . i_ 

9t ij 9t 3 ij 

life can show that 


JL^R -_J i_R 


(4.13.2) 


9 

9t 




D - 2 -ii (^+ $ 

^ pq 



(4.13.3) 


where is the production tensor, 4>jLj ns the pressure-strain 

tensor, and as the dissipation tensor. Ihe turbulent kinetic 

energy production is (P , the turbulent dissipation is e, and the 


80 



trace of the pressure strain tensor is $ . The trace of the Reynolds 

2 

stress tensor is pq . 

Setting (8/3t)b^j equal to zero and solving for j , we have 




ij 


- ‘■ij> - - “ij> 


+ 2 $ 






(4.13.4) 


where we have used for the production anisotropy, defined as 


Ij 


3 ij 


(4.13.5) 


Equation (4.13.4) gives us a form from which to construct models. To 
test the assumption that (3/3t)bj^j = 0, we calculated the ratio of the 
right side of (4.13.4) to the left side and calculated the average value 
and variance of this ratio. If (4.13.4) were identically satisfied, 
this ratio would be one and the variance would be zero. Table 4.14 
shows these results . Wfe see that this assumption is reasonable for 
the i = 1, j = 1 and i = 3, j = 3 components, but questionable for 
the other two. 

life constructed several models by simplifying. Ignoring the trace 
of the pressure-strain tensor and using the observed fact that the dis- 
sipation anisotropy is proportional to the stress anisotropy, we have 


= 2(P(h^. - p^.) - 2ebjj(l-c^) (4.13.6) 

where Cj represents this proportionality. Ihe average value of Cj^ 
is shown in Thble 4.15, along with its variance. Ife also show the re- 
sults graphically in Fig. 4.29, in the same way we presented the values 
in Section 4.9. 

There is reasonable agreement for the diagonal terms, but the off- 
diagonal or shear-stress term is not well represented. Since cj^ rep- 
resents the proportionality between b^^j and we know that its 
value should be about one . Ife suspect that the failure on the off- 
diagonal term is due to the changing value of ^ seen ear- 
lier, this off-diagonal compnent of the stress anisotropy tensor changes 
more rapidly throughout the "good time" than the diagonal components . 
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The failure in this term may be due to the lack of structural similarity 
in the shear stress . 

Searching for improvements to this model, we put an adjustable con- 
stant in front of the first term in (4.13.6) and rewrite the equation as 

*^ij “ ^l^^^ij “ Pij^ “ (4.13.7) 

where the second term is written to appear like the Rotta model. Note 
that (4.13.7) looks very much like the Gibson-Launder plus Rotta models, 
but with the Inclusion of in the first term. Table 4.16 shows the 

values of c^^ and C 2 determined from a least-squares fit using 

(4.13.7) as the fitting function, fcfe again see a large variation in the 
values for 1=3, j = 3. Ihe agreement among the other three equations 
is not good, with a four-to-one variation on C 2 and two-to-one on Cj^ . 

Perhaps models based on structural similarity would be appropriate 
for this flow, if we could carry the simulation further in time. How- 
ever, we are limited to simulation times before the flow reaches this 
hypothetical similarity, and therefore we should not dismiss this idea. 
It deserves a reinvestigation in a more fully developed homogeneous 
shear flow. 
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Table 4.1 


NON-DIMENSIONAL GROUPS 



Computational Range 

Experimental Range 

SL 



q 

CT 

8.5 < — < 26.6 

7.8 < -^ < 14.8 

Shear Number 

q 

q 

M = 1 

c 

Fluctuating Mach 
Number 

0.06 < M < 0.31 

< 0.3 

Re, 



A p 

18.4 < Re;^ < 120.6 

133 < Re;^ < 398 

Taylor Microscale 
Reynolds Number 




Table 4.2 

TABLE OF SIMULATIONS 



D«jsignation 

Type 

Parameter Range 

SL/q 

M = q/c 



IH64A 

Isotropic 

0.0 

0.078-0.034 

15.0-40.0 

A. 

HS64A 

Shear 

22.0-25.8 

0.064-0.065 

18.4-23.3 

B. 

HS64B 

Shear 

9.7-12.5 

0.144-0.146 

30.1-39.0 

C. 

HS64C 

Shear 

8.5-11.1 

0.312-0.316 

40.7-53.3 

D. 

HS64D 

Shear 

13.1-14.9 

0.207-0.222 

43.6-56.7 


HS64F 

Shear 

24.0-26.6 

0.238-0.262 

54.9-77 .7 


HS64G 

Shear 

9.9-13.3 

0.268-0.273 

29.3-37.9 

H. 

HS64H 

Shear 

10.7-11.4 

0.250-0.282 

93.8-120.6 
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Table 4.3 


INVARIANTS OF THE REYNOLDS STRESS ANISOTROPY TENSOR b^j 
Fitting function: f = d (1+bM^) (Re^) ^ 


II : 


a 

b 

c 

d 


0.54 

0.61 

0.19 

0.01 


rms error = 0.048 


HI : 


a 

b 

c 

d 


1.27 

2.37 

0.55 

0.001 


rms error = 0.010 


Table 4 .4 

INVARIANTS OF THE DISSIPATION ANISOTROPY TENSOR d 

/ CT 0 

Fitting function: f = d (1+bM )(Rej^) 


ij 


II : 


a 

b 

c 

d 


0.688 
-2 .079 
0.334 
0.008 


err = .072 


III : 


a 

b 

c 

d 


1.787 

-6.142 

1.473 

0.001 


err = 0.024 
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Table 4.5 


ROTTA MODEL EVALUATION 

C T ^ O A 

Fitting function: f = d (1+bM )(Rfij^) 


Rotta Tern Alone (4’ij) 


Equation 

Cl 

a 

b 

c 

d 

rms 

Error 

1 = 1, :i = 1 

1.114 * 0.723 

-1.519 

-3.763 


1.949 

0.049 

1=2, :i = 2 

1.730 ± 0.899 

-0.776 

3.614 


0.762 

0.074 

1 = 3, j = 3 

0.680 ± 0.583 

-1.965 

-7.352 

1.226 

1.481 

0.119 

1 = 1, :i = 2 

1.559 ± 1.090 

-0.592 

-2.311 

1.219 

0.067 

0.079 

Rotta Term with 

Dissipation Anisotropy 

(<l>ij-2ed 




1 = 1, J = 1 

2.886 ± 0.722 

-0.458 

-2 .840 


2.059 

0.057 

1 = 2, j = 2 

3.727 ± 0.824 

-0.352 

0.105 

0.325 

2.652 

0.049 

1 = 3, ;i = 3 

1.719 ± 0.921 

-0.126 

-4.981 

1.059 

0.053 

0.154 

1 = 1, j = 2 

3.035 ± 1.122 

-0.294 

-2.252 

0.725 

0.445 



0.057 


Table 4.6 

CONSTANTS IN LUMLEY'S FITTING FUNCTION 
FOR THE ROTTA CONSTANT 
(See Eq. (4.9.6)) 


Equation 

a 

b 

c 

d 

rms 

Error .. 

BIB 

7.795 

244.196 

21.233 

-5.81 

0.067 


4.275 

152.380 

21.765 

-1.842 

0.043 


9*191 

525.294 

2.106 

31 .064 

0,054 


Lumley predicts a = 7.77, b = 80.1, c = 62.4, d = 2.3. 
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Table 4.7 


LUMLEY AND NEMAN'S FITTING FUNCTION 
FOR THE ROTTA CONSTANT 

f = 2.0 + a Re^'^^ 

(See Eq. (4.9.8)) 


Equation 

a 

rms Error 

i = 1, j = 1 

1.682 

0.334 

i = 2, j = 2 

3.681 

0.376 

i = 3, j = 3 

-1.550 

0.411 

. i = 1 , 3 = 2 

1.624 

0.439 


Table 4 .8 

GENERAL FOURTI^RANK TENSOR MODEL 
Fitting function; f = d (“) (1+bM^) (Re^) 

(See Eq. (4.9.8)) 


Equation 

-A 

a 

b 

c 

d 

rms 

Error 

i = 1, j = 1 

-4.090 ± 0.173 

0.088 

-0.225 

0.061 


0.010 

1 = 2, j = 2 

-2.134 ± 0.052 

0.073 

0.799 

-0.001 


0.014 

i = 3, j = 3 

-0.482 db 0.147 

-0.701 

-3.593 

-0.333 

-13.190 

0.088 

1 1. J = 2 : 

—1, 153 ± 0.098 

0.190 

-0.675 

0.030 

-0.642 

0.025 
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Table 4.9 


GIBSON /LAUNDER MODEL FOR FAST TERM 
Fitting function: f = d (l+bM^)(Re^)^ 


Equation 

^2 

a 

b 

c 

d 

rms 

Error 

i = 1, j = 1 

0.324 ± 0.061 

-0.423 

0.707 

-0.274 

2.673 

0.047 

i = 2, j = 2 

-0.134 ± 0.052 

0.471 

48.117 

-0.156 

-0.020 

0.179 

1=3, j = 3 

0.289 ± 0.088 

-0.701 

-3.592 

-0.333 

7.914 

0.088 

1=1. j = 2 

0.198 ± 0.029 

0.000 

-0.705 

-0.236 

0.501 

0.059 


Table 4.10 

GENERAL FOURTH-RANK TENSOR WITH 
DISSIPATION ANISOTROPY TERMS FOR FAST TERM 


Equation 


C2 

C3 

rms 

Error 

i = 1, j = 1 

-1.818 

1.624 

0.612 

0.033 

i = 2, j = 2 

-2.221 

1.624 

0.612 

0.016 

i = 3, j = 3 

-1.926 

0.403 

0.0 

0.217 

i = 1, j = 2 

-1.731 

0.132 

-0.055 

0.011 
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Table 4.11 


COMPRESSIBLE PRESSURE-STRAIN TERMS 
General Fourth-Rank Tensor Model 
(Compressible terms alone) 

Si 

Fitting function: f = d| m'^(Roj^)^ 


Equation 

^1 

a 

b 

c 

d 

rms 

Error 

i = 1, j = 1 

-3.962 ± 0.199 

-0.015 

-0.111 

0.095 

-2.397 


i = 2, j = 2 

^1.989 ± 0.076 

-0.101 

-0.083 

0.034 

-1.995 


i = 3, j = 3 

-0.678 ± 0.138 

-0.430 

0.071 

-0.232 

-5.663 

0.093 

i =1, j = 2 

-1.228 ±0.094 

0*208 

0*060 

-0*033 

-0.882 

0.033 


Table 4.12 

COMPRESSIBLE PRESSURE-STRAIN TERMS 
General Fourth-Rank Tensor Model on Compressible and 
Fast Pressure-Strain Terms 

Fitting function: f = d (1+bM^) (Re^)^ 


Equation 

h 

a 

b 

c 

d 

rms 

Error 

1 = 1. j = 1 

-3.052 ± 0.323 

0*121 

-1.889 

0.167 



1=2, j = 2 

-2.123 ± 0.068 

-0.054 

-0.065 

0.009 



i = 3, j = 3 

-1.160 ± 0.261 

-0.556 

-1.232 

-0.262 

-14.409 



-1*388 ± 0*069 . 

0^073. 

.0..592 

-0*065. 

^1.418- 
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Table 4.13 


SIMULTANEOUSLY DETERMINED CONSTANTS FOR THE TOTAL 
PRESSURE-STRAIN TENSOR PLUS DISSIPATION ANISOTROPY 

Modeling of ~ 

Gibson/Launder + Rotta Models 


Equation 

Fast Const 

Rotta const 

rnis Error 

i = 1, j = 1 

0.334 

3.234 

0.120 

i = 2, j = 2 

0.022 

3.106 

0.099 

i = 3, j = 3 

1.122 

-4.571 

0.163 

i = 1, j = 2 

0*426 

0.899 

- ■ - 

0.057 


Fourth-Rank Tensor + Rotta Models 


Equation 

Fast const 

Rotta const 

rms Error 

i = 1, j = 1 

-3.332 

3.234 

0.325 

i = 2, j = 2 

-1.978 

3.106 

0.042 

i = 3, j = 3 

-1.869 

-4.571 

0.163 

i = 1, j = 2 

-1.500 

2.074 

0.020 
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Table 4.14 


TEST OF STRUCTURAL SIMILARITY 


"l = 




-ij) - - “ij) + *C=ij ^ V 




Fitting function: f = d (1+bM^) (Re^^)*^ 


Equation 

Cl 

a 

b 

c 

d 

rms 

Error 

i = 1, j = 1 

0.868 =t 0.067 

-0.120 

0.475 

0.045 

0.979 

0.045 

CM 

11 

CM 

II 

•H 

0.616 ± 0.208 

-0.078 

10.25 

-0.062 

0.618 


i = 3, j = 3 

0.944 ± 0.09 

-0.207 

-0.862 

0.089 

1.214 


1=1. j = 2 

1.381 ± 0.268 

0.379 

0 .403 

0.212 

0.218 

0.099 


Table 4.15 

STRUCTURAL SIMILARITY MODEL 
“^^Ij ^ “ 2eb^j(l-Cj^), see Eq . (4.12.5) 


Fitting function: f = d (~) (1+bM^) (Re^^)^ 


Equation 

Cl 

a 

b 

c 

d 

rms 

Error 

■i 

1.005 ± 0.170 

-0.012 

-2.811 

-0.038 

Mj 

0.125 

nM 

1.340 ± 0.202 

-0 .069 

-2.564 

0.218 


0.118 

|HB 

0.920 ± 0.452 

0.873 

3.833 

0.006 


0.318 

■HB 

0.023 ± 0.529 

10.570 

-224.900 

-4.53 

0.0 

0.649 
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Table 4.16 


STRUCTURAL SIMILARITY MODEL 


Equation 


C2 

rms Error 

i = 1, j = 1 

0.977 

0.128 

0.057 

1 = 2, j = 2 


0.217 


i = 3, j ■= 3 


-0.946 


: 1 ■ = , 1 , , j = 2 ■ , , 

1.454 

0,245 

0.204 
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Chapter V 

CONCLUSIONS AND RECOMMENDATIONS 

While much turbulence research has been done In incompressible 
flows, most turbulent flows of technological interest are compressible. 
Knowledge about turbulence modeling gained from incompressible studies 
has been widely used in compressible simulations, sometimes with excel- 
lent results and other times with failure . Because measurements in 
high-spe( 2 d flows are so difficult, little is known about the structure 
of the Rfjynolds stresses and their dynamic equations in these cases, and 
consequently there is little guidance on how to construct turbulence 
models . 

To study these stresses , we used the power of a large , modern 
vector computer to simulate directly the full, compressible Navier- 
Stokes equations, with no turbulence model. By using the computer in 
this way, as a "numerical wind tunnel," we were able to measure turbu- 
lence quantities that are of particular interest to modelers . Some of 
these quantities, for example the pressure-strain terms in the Reynolds 
stress equations, cannot be directly measured by experimental means. 
Tltiis type of simulation allows these terms to be directly measured and 
studied for the first time. 

To this end, we chose a simple, homogeneous shear flow that is an 
approximation to a small part of a more complicated shear flow. Wb have 
extended and developed incompressible techniques for the simulation of 
compressible homogeneous flows with a general mean-velocity gradient. A 
code that runs on the ILLIAC IV computer was constructed to implement 
these techniques. After thorough testing, eight complete 64 x 64 x 64 
mesh size simulations were performed, and the resulting simulated flow 
fields were used as a data base in which to study turbulence quantities . 

Several measures of the structure of the Reynolds stresses were 
used to compare the simulated stress tensor to measurements from experi- 
ment. Agreement was good, with little indication that the stresses in 
compressible shear flows are significantly different from the incompres- 
sible case . 
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A spectral decomposition originally due to Moyal (1951) was im- 
plemented and extended numerically to divide the velocity field into 
"incompressible" and "compressible" parts. Vfe found that the stresses 
arise primarily from the incompressible part of the velocity field, in 
agreement with Ifarkovin's (1962) hypothesis that compressibility has but 
a small effect on these stresses. 

Although the question of whether or not a homogeneous shear flow 
comes to a state of structural equilibrium is yet unresolved, we find an 
indication that these shear flows are tending toward similarity, as 
defined by constant values of the stress anisotropy tensor. Vb are 
unable to carry the simulations far enough in time to see this state, 
but the time histories of the anisotropies seem to indicate this. This 
has an important bearing on the construction of some of the turbulence 
models that we have evaluated . 

It is known that the performance of incompressible turbulence 
models deteriorates at high Mach number. Ib search further for this 
effect, we studied the Reynolds stress equations and evaluated models 
for the various terms . 

In a compressible flow, the trace of the pressure strain tensor is 
not zero. M.thout direct measurement, it cannot be told from the 
turbulent kinetic energy equation whether this term is productive or 
dissipative. Ivb find it to be dissipative. This is an effect that is 
not represented in current kinetic energy equation models , but perhaps 
it should be . 

To gain insight into pressure strain-term modeling in incompress- 
ible flows, these terms are symbolically decomposed into two parts, the 
Rotta and the Fast terms. This decomposition is not actually performed 
but is simply used to Justify the construction of models for each part . 
Ife performed this decomposition numerically for the first time and were 
able to evaluate the individual terms directly. 

In compressible flow we found that this decomposition produces 
three parts. Because the third part represents the deviation of the 
pressure-strain tensor from the incompressible case, we call this third 
part the compressible term. It was found to be of the same order of 
magnitude as the Rotta and Fast terms, and therefore it must be modeled. 
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Models were separately evaluated for each part of the pressure 
strain tensor. Vb found that the Rotta model performs much better when 
used to replace the sum of the Rotta pressure-strain term and the dis- 
sipation anisotropy, and would recommend that it be used this way. 

Several models were evaluated for the Fast and compressible terms, 
with varying degrees of success. Although its performance must be 
significantly improved, the general linear fourth-rank tensor model 
performs the best. There are significant mathematical objections to its 
use, but until more sophisticated models are developed we would recom- 
mend its use, with carefully chosen constants for the sum of the Fast 
and Compressible terms . 

Ife found that the constants that must be determined for use in 
these models are more uniform if they are determined simultaneously from 
the total pressure strain tensor. Wfe suspect that individual determina- 
tion of the constants from the decomposed pressure-strain field is not 
advisable and recommend their simultaneous determination . 

A class of models based on structural similarity was proposed for 
the entire pressure-strain tensor. These models show promise but need 
to be studied in a more fully developed turbulent flow. 

The dissipation tensor was found to be very anisotropic . I* 
searched for a dependence of this anisotropy on the Reynolds number and 
found a positive relation between the two, indicating that the dissipa- 
tion is not necessarily isotropic at higher Reynolds number in a per- 
sistently sheared flow, as was thought. This important result may alter 
thinking about the structure and the modeling of the dissipation terms, 
vrhich are normally considered to be isotropic at high Reynolds numbers . 

The results of this work have fallen primarily into two categories : 
results that bear on all turbulent shear flows and results that pertain 
only to compressible shear flows. Vb have identified some differences 
that occur in compressible flows and evaluated models for both new and 
old terms. However, much that we have studied applies to all shear 
flows . 

Vfe have concentrated primarily on the Reynolds stresses, but there 
is much additional work that needs to be done. 


94 



For further research In this area, we recommend that: 

1. The turbulent terms in the thermal energy equation be studied and 
models compared with the exact flow fields . 

2. The capability be developed to follow the flow for a greater non- 
dimensional time (a larger computer?). 

3. An exactly parallel Incompressible code be developed to better 
assess the changes in the compressible flow. 

4. The method be extended for simulation of higher Reynolds number 
flows (possibly via large eddy simulations?). 

5. The method be extended for non-periodic boundary conditions (a 
compressible mixing layer would be very interesting for noise 
studies .) 

This work simply scratches the surface of the vast new area of 
direct simulation of turbulent flows. Use of the computer as a numer- 
ical wind tunnel complements the use of laboratory wind tunnels and 

allows us to study turbulence in a way that has never before been pos- 
sible . It will undoubtedly give us great insight when the next genera- 
tion of super-computers becomes available, as it already does right now. 
It is the author's opinion that this technique will become a major tool 
in the turbulence modeler's inventory in the next few years. 
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MEAN FLOW 



1.1. Schematic of the Mean Velocity Profile. 




Fig. 2.1. Modified wave numbers for two finite- 

difference approximations and the pseudo- 
spectral approximation to spatial deriva- 
tives vs. the analytic wave number. 


Fig. 2.2. Computational root of the Euler explicit 
method for linear Burgers equation in the 
complex plane. Viscous stability number 
V = 0 (convective terms only) , □ and the 
analytic root, Q ; for several Courant 
numbers . 
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Fig. 2.3. 


Computational root of the Euler explicit Fig. 2.4. 
method for linear Burgers equation, Courant 
number C = 0 (viscous terms only) , and 
the analytic root vs. the viscous stability 
number V = hk^/v. 


Argument (phase angle) of the computational 
root of the Crank-Nicholson method vs. the 
argument of the analytic root. Linear Bur- 
gers equation; viscous stability number 
V = 0 (convective terms only) . 
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1.0 



Fig. 2.5. Roots for analytic solution and Crank- 
Nicholson method vs. viscous stability 
number. Linear Burgers equation, Courant 
number C = 0 (viscous terms only) . 


Fig. 2.6. Computational root of MacCormack's method in 
the complex plane. Viscous stability number 
V = 0 (convective terms only) . For several 
wave numbers. Note the highly dissipative 
( [aj << 1.0) behavior for the higher wave 
numbers. The wave numbers are labeled kA =, 
the Courant numbers are labeled C =. The 
analytic root Q liss on the unit circle and 
is labeled with its Courant number. 
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Fig. 2. 



7. Roots for analytic solution and 4th-order 
Runge-Kutta method in the complex plane. 
Linear Burgers equation. Viscous stability 
number V = 0 (convective terms only) , for 
several Courant numbers. Analytic root 0> 
Runge-Kutta root 0 . 


Fig. 2.8. Roots for analytic solution and 4th-order 

Runge-Kutta method. Linear Burgers equation. 
Courant number C = 0 (viscous terms only) . 
vs. viscous stability number V. Analytic 
root O, Runge-Kutta root Q. 



64 pts. 



Fig. 3.1. The pencil system data base method for handling 
movement of data between the ILLIAC IV disk and 
core memory. Data are brought into core as a 
"pencil" which is an 8 x 16 x 64 point section 
of the 64 X 64 X 64 size mesh. Each mesh point 
is part of three pencils corresponding to the 
X, y, and z directions. 
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Fig. 3.2. The ''remeshing" process for the computational 
mesh. Top figure: St = 0.0 (mesh metric) Car- 
tesian. Middle figure: St = 0.5. Computation 
is stopped here and the flow field is interp- 
olated onto the mesh in the bottom figure at 
St = -0.5. Computation is then resumed. 
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3-D ENERGY SPECTRA 



Fig. 3.3. 


3-D spectra of the Initial velocity field 
for simulation IH64A. All Initial spectra 
are similar except for the intensity of 
the turbulence. 

E indicates twice the turbulent kinetic 


El 

E2 

E3 


energy pu^u^. 
indicates pu^. 
indicates pv^. 
indicates pw^. 


Wave No. indicates the magnitude of the 


wave number k. 
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3-D ENERGY SPECTRA 
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El 

E2 

E3 



WAVE MO. 


Fig. 3.4. 3-D spectra of the velocity Meld for the 

isotropic homogeneous simulation IH64 at 
time = 7.849. Same nomenclature as for 
Fig. 3.3. 







I wo HOiNi OUKKLLAliONS 


1— D ENERGY SPECTRA 



Fig. 3.5. Two-point correlations in the z-direction Fig. 3.6. 1-D spectra of tlie velocity field for simu- 

for simulation IH64A at time = 7.849. lation IH64A at time = 7.849. Uave No. 

Rll (z) indicates the u velocity corr. indicates the wavenumber in the re- 

R22 (z) indicates the v velocity corr. spective direction. 

R33 (z) indicates the w velocity corr. 

Delta R indicates the separation (Delta R 

=3.14 is half the computational box size.) 
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MICRO SCALES 


INTEGRAL SCALES 



Fig. 3.7. z-direction Taylor microscales vs. time 
for simulation IH64A. 


U(Z) 

indicates 

the 

scale 

assoc, with 

u 

V(z) 

indicates 

the 

scale 

assoc, with 

V 

W(Z) 

indicates 

the 

scale 

assoc, with 

w 


Fig. 3.8. z-direction integral length scales vs. 

time for simulation IH64A. Same nomen- 
clature as for Fig. 3.7. 




ENERGY HISTORY 


SKEWNESS HISTORY 



TIME time 


Fig. 3.9. Twice the turbulent kinetic energy and Fig. 3.10. Velocity derivative skewnesses vs. time 

turbulent intensities vs. time for slmu- for simulation IH64A. 

lation IH64A. SKl(K) is the skewness assoc, with Uj^ . 

Same nomenclature as for Fig. 3.3. SK2(K) is the skewness assoc, with u^2,2' 

SK3 (K) is the skeT>mess assoc . with u^ ^ • 
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STRESS COEFFICIENTS 


ENERGY HISTORY 



S-TIME S*TIME 


Fig. 3.11. Shear stress correlation coefficients vs. Fig. 3.12. Twice the turbulent kinetic energy and the 

St for simulation HS64B. The uv corre- turbulent Intensities vs. time for the 

lation is defined in Eq. (3.7.1). simulation HS64B. Same nomenclature as 

for Fig. 3.3. 




3-D ENERGY SPECTRA 


3-D ENERGY SPECTRA 









1-D ENERGY SPECTRA 




WAVE NO. 


Fig. 3.15. 1-D spectra of the velocity field for 

HS64B at St = 6.0. WAVE NO. indicates 
the wave number in the respective 

directions . 


TWO POINT CORRELATIONS 



DELTA R 


Fig. 3.16. Two-point correlations in the z-direction 
for HS64B at St = 4.0. Nomenclature same 
as for Fig. 3.5. 
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TWO POINT CORRELATIONS IWU PUiNi CukRELATIONS 



DELTA R delta R 


Fig. 3.17. Two-point correlations in the z-direction Fig. 3.18. Two-point correlations in the z-direction 
for HS64B at St = 6.0. Nomenclature for HS64B at St = 7.0. Nomenclature 

same as for Fig. 3.5. same as for Fig. 3.5. 
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TWO POINT CORRELATIONS 


PARAMETER SPACE SL(11)/Q VS. AVRG. MACH NUM. 



Fig. 4.1. Definitions of integral length scales 

illustrated for a typical two-point cor- 
relation. The arrow at 0.42 indicates 
Corrsin's definition. The arrow at 1.42 
indicates our definition. 


Fig. 4.2. Parameter space of shear number SL(11)/Q 
vs. fluctuating Mach number for the seven 
shear-flow simulations. 
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Fig. 4.3. Parameter space of shear number SL(11)/Q 
vs. Taylor microscale Reynolds number for 
the seven shear-flow simulations. 


Shear stress correlation coefficients 
vs. St. Letters in the figure number 
correspond to the individual shear-flow 
s imulat ions . 


Fig. 4.4a 
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STRESS COEFFICIENTS 


STRESS COEFFICIENTS 




Fig. 4.4.b. Shear stress correlation coefficients 
vs. St. Letters in the figure number 
correspond to the individual shear- 
flow simulations. 


Fig. 4.4 b. Shear stress correlation coefficients 
vs. St. Letters in the figure number 
correspond to the individual shear-flow 
simulations . 
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STRESS COEFFICIENTS 


STRESS COEFFICIENTS 



Fig. 4.4d. Shear stress correlation coefficients 
vs. St. Letters in the figure number 
correspond to the individual shear- 
flow simulations. 


Fig. 4.4f. Shear stress correlation coefficients 
vs. St. Letters in the figure number 
correspond to the individual shear-flow 
simulations . 
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STRESS COEFFICIENTS 


STRESS COEFFICIENTS 



Fig. 4.4g. Shear stress correlation coefficients 
vs. St. Letters in the figure number 
correspond to the individual shear-flow 
simulations. 


Fig. 4.4h. Shear stress correlation coefficients 
vs. St. Letters in the figure number 
correspond to the individual shear flow 
simulations . 
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Fig. 4.5. Angle, a, between Cartesian and princi- Fig. 4.6. 
pal axis coordinates of the Reynolds stress 
tensor vs. estimated values of a from the 
least -squares fitted function. The HGC 
flow is indicated by the solid dot. 


Ratio of the principal stresses of the Rey- 
nolds stress tensor in the x-y plane vs. 
estimated values from the least— squares 
fitted function. The HGC flow is indicated 
by the solid dot. 





REYNOLDS STRESS ANISOTROPY REYNOLDS STRESS ANISOTROPY 



S*TIME . S'TIME 


Fig. 4.7a, Reynolds stress anisotropy, ^^8- 4.7b. 

St for the seven shear simulations 
bi2 = RUV, b^^ = RUU, b22 = RVV, b23 = 

RWW. Letters in the figure numbers cor- 
respond to the individual shear-flow 
simulations . 


Reynolds stress anisotropy, vs. 

St for the seven shear simulations 

hi " ^11 " ^22 " ^33 " 

RWW. Letters in the figure numbers cor- 
respond to the individual shear-flow 
simulations . 




- 0.25 - 0.20 



S*TIME S*TIME 


Fig. 4.7c. Reynolds stress anisotropy, , vs. Fig. 4.7d. Reynolds stress anisotropy, vs. 

St for the seven shear simulations, St for the seven shear simulations 

bi2 = RUV, = RUU, b^2 = RW, b^^ = ^12 " ^22 ^33 " 

RWW. Letters in the figure numbers cor- RWW. Letters in the figure numbers cor- 
respond to the individual shear -flow respond to the individual shear-flow 

simulations. simulations. 







REYNOLDS STRESS ANISOTROPY 


REYNOLDS STRESS ANISOTROPY 



Fig. 4.7f. Reynolds stress anisotropy, 

St for the seven shear simulations 

= RUV, = RUU, b^2 = RVV, b33 = 

RWW. Letters in the figure numbers cor- 
respond to the individual shear -flow 
simulations . 


Fig. 4.7g. Reynolds stress anisotropy, ^ij > vs. 

St for the seven shear simulations 

^11 " ^22 " ^33 " 

RWW. Letters in the figure numbers cor- 
respond to the individual shear-flow 
simulations . 
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RLYNOLDS SiKLSb AN i so I ROPY 


LUMLEYS INVARIANT MAP 




Fig. 4.7h. Reynolds stress anisotropy, 

St for the seven shear simulations. 
b^2 = RUV, = RUU, b^2 = RVV, b33 

= RWW. Letters in the figure numbers 
correspond to the individual shear-flow 
simulations . 


Fig. 4.8. Invariants of the Reynolds stress anisotropy. 
II vs. Ill (defined in Eq. (4.5.6). All 
turbulence must be in the triangular region. 
The left and right lines indicate axisymmetric 
turbulence. The top line indicates two- 
dimensional turbulence. The HGC flow is 
indicated by the black dot. 
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DECOMPOSED RUU AND RWW 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 00 0 5 1-0 1.5 2.0 2.5 3.0 35 4.0 4.5 5 0 5.5 6.0 

S*TIME S*TIME 


Fig. 4.9a. Moyal decomposed Reynolds stresses Fig. 4.9b. Moyal decomposed Reynolds stresses R^^^ 

and R-., as defined by Eq. (4.6.3) and R-., as defined by Eq. (4.6.3) 

vs. Str vs. Stf-" 

RUU SOL = R^j^, RUU DIV = R°^, RWW SOL = R^^. RWW DIV = R°^ 

Letters in the figure numbers correspond to the individual shear flow 
simulations. 





DECOMPOSED RUU AND RWW 
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S*TIME 
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Moyal decomposed Reynolds stresses R-, Fig. 4.9d. Moyal decomposed Re 3 molds stresses 
and R,_ as defined by Eq. (4.6,3) and defined by Eq. (4.6.3) 

vs. StP vs. St. 

RUU SOL = R^.,, RUU DIV = RWW SOL = RWW DIV = R^^ 

Letters in the figure numbers correspond to the individual shear flow 
simulations . 





Fig. 4.9f. Moyal decomposed Re3molds stresses Fig. 4.9g. Moyal decomposed Reynolds stresses R^^ 

and R„. as defined by Eq. (4.6.3) and R„„ as defined by Eq. (4.6.3) 

vs. St. vs. St. 

RUU SOL = R^^, RUU DIV = R^^^, RWW SOL = RWW DIV = R^3 

Letters in tbe figure numbers correspond to the individual shear flow 
simulations. 
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decomposed RUU and RwW 



Fig. 4.9h. Moyal decomposed Reynolds stresses R^^ 
and R ^2 as defined by Eq. (4.6.3) 
vs. St. „ Q 

RUU SOL = RJ^ , RUU DIV = R, , , RWW SOL = 

R^^, RVJl-J DIV = Rg^" Letters in the fig- 
ure numbers correspond to the individual 
shear-flow simulations. 


Fig. 4.10a. Moyal decomposed Reynolds stresses R^^^ 
and R|^^ defined by Eq. (4.6.3) 

vs. St. « TV 

RUV SOL = R, , RUV DIV = IL , RVV SOL ^ 

^ 22 ’ ~ ^ 22 ’ 

ure numbers correspond to the individual 
shear-flow simulations. 
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DECOMPOSED RUV AND RVV 


DECOMPOSED RUV AND RVV 



S.TIME 1,. S-TIME 


S D S' 

Fig. 4.10b. Moyal decomposed Reynolds stresses R^’ Fig. 4.10c. Moyal decomposed Reynolds stresses R^’ 

and ^22^ defined by Eq. (4.6.3) and R 22 defined by Eq. (4.6.3) 

vs. St. vs. St. 

RUV SOL = " ^2’ " ^ 22 ’ " ^2 

Letters in the figure numbers correspond to the individual 
shear-flow simulations. 
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DECOMPOSED RUV AND RVV 



Fig. 4.10d, Moyal decomposed Reynolds stresses R^’° Fig. 4.10f. Moyal decomposed Reynolds stresses R^^’ 
and R 22 ° as defined by Eq. (4.6.3) - and r|’^ as defined by Eq. (4.6.3) 

vs. St. vs. St. 

RUV SOL = R^2’ " ^2’ " ^22’ " ^2 

Letters in the figure numbers correspond to the individual shear- 
flow simulations. 
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Fig. 4.10g. Moyal decomposed Re3niolds stresses Fig- 4.10h. Moyal decomposed Reynolds stresses 

^12^ and ^22^ defined by Eq. ^22^ defined by Eq. 

(4.6.3) vs. St. (4.6.3) vs. St. 

RUV SOL = R^2» = ^2’ ^ 22 * 

Letters in the figure numbers correspond to the individual shear-flow simulations. 
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SOLLNOIDAL LNLKGY iPLCTRA 


SOLENOIDAL ENERGY SPECTRA 




Fig. 4.11a, 3-D spectra of the solenoidal part of Fig. 4.11b. 3-D spectra of the solenoidal part of 
the Reynolds stress field as defined the Reynolds stress field as defined 

by Eq. (4 .6.5) . by Eq. (4.6.5) . 

Nomenclature same as for Fig. 3.3. Letters in figure numbers correspond to the individual 
shear-flow simulations . 
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SOLENOIDAL ENERGY SPECTRA 


SOLENOIDAL ENERGY SPECTRA 




Fig. 4.11c. 3-D spectra of the solenoidal part of Fig. 4. lid. 3-D spectra of the solenoidal part of 
the Reynolds stress field as defined the Reynolds stress field as defined 

by Eq. (4.6.5). by Eq. (4.6.5). 

Nomenclature same as for Fig. 3.3. Letters in figure numbers correspond to the individual 
shear-flow simulations. 
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Fig. 4. Ilf. 3-D spectra of the solenoidal part of Fig. 4.11g. 3-D spectra of the solenoidal part of 
the Reynolds stress field as defined the Reynolds stress field as defined 

by Eq. (4.6.5). byEq. (4.6.5). 

Nomenclature same as for Fig. 3.3. Letters in figure numbers correspond to the individual 
shear-flow simulations. 
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SOLENOIDAL ENERGY SPECTRA DIVERGENCE ENERGY SPECTRA 



WAVE NO. WAVE NO. 


Fig, 4.11h. 3-D spectra of the solenoidal part of Fig. 4.12a. 3-D spectra of the dilatation part of the 
the Reynolds stress field as defined Reynolds stress field as defined by Eq. 

by Eq. (4.6.5). Nomenclature same as (4.6.5). Nomenclature same as for Fig. 

for Fig. 3.3. Letters in figure numbers 3.3. Letters in figure numbers correspond 

correspond to the individual shear-flow to the individual shear-flow simulations, 

simulations. 
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DIVERGLNCL lNlkGY spectra 


DIVERGENCE ENERGY SPECTRA 




Fig. 4. 12b. 3-D spectra of the dilatation part of Fig. 4.12c. 3-D spectra of the dilatation part of the 

the Reynolds stress field as defined Reynolds stress field as defined by Eq. 

by Eq. (4.6.5). (4.6.5) 

Nomenclature same as for Fig. 3.3. Letters in figure numbers correspond to the individual 
shear-flow simulations. 
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DIVERGENCE ENERGY SPECTRA 


DIVERGENCE ENERGY SPECTRA 




Fig. 4.12d. 3-D spectra of the dilatation part of Fig. 4.12f. 3-D spectra of the dilatation part of 
the Reynolds stress field as defined the Reynolds stress field as defined 

by Eq. (4.6.5). by Eq. (4.6.5). 

Nomenclature same as for Fig. 3.3. Letters in figure numbers correspond to the individual 
shear-flow simulations. 
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DIVERGENCE ENERGY SPECTRA 


DIVERGENCE ENERGY SPECTRA 




Fig. 4.12g. 3-D spectra of the dilatation part of Fig. 4.12h. 3-D spectra of the dilatation part of 
the Reynolds stress field as defined the Reynolds stress field as defined 

by Eq. (4.6.5). by Eq. (4.6.5). 

Nomenclature same as for Fig. 3.3. Letters in figure numbers correspond to the individual 
shear -flow simulations. 
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DECOMPOSED STRESS SPECTRA 


DECOMPOSED STRESS SPECTRA 



Fig. 4.13a. 3-D spectra of the Moyal decomposed Fig. 4.13b. 3-D spectra of the Moyal decomposed 

shear stress as defined by Eq. (4.6.5). shear stress as defined by Eq. (4.6.5). 

S D 

V1V2(K) indicates Rj^ 2 » V1W2(K), W1V2(K), WiW2(K) are the contributions to 

WAVE NO. indicates the magnitude of the wavenumber vector. Letters in the figure 
numbers correspond to the individual shear flow simulations. 







DECOMPOSED STRESS SPECTRA 


DECOMPOSED STRESS SPECTRA 



Fig. 4.13c. 3-D spectra of the Moyal decomposed Fig. 4.13d. 3-D spectra of the Moyal decomposed 

shear stress as defined by Eq. (4.6.5). shear stress as defined by Eq. (4.6.5). 

V1V2(K) indicates V1W2(K) , W1V2(K), WiW2(K) are the contributions to R^^* 

WAVE NO. indicates the magnitude of the wavenumber vector. Letters in the figure 
nimibers correspond to the individual shear flow simulations. 
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Fig. 4.13f. 3-D spectra of the Moyal decomposed Fig. 4.13g. 3-D spectra of the Moyal decomposed 

shear stress as defined by Eq. (4.6.5). shear stress as defined by Eq. (4.6.5) 

V1V2(K) indicates ^^2’ V1W2(K), W1V2(K), W1W2(K) are the contributions to ^ 2 ' 

WAVE NO. indicates the magnitude of the wavenumber vector. L etters in the figure 
ntambers correspond to the individual shear-flow simulations. 
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Fig. 4.13h. 3-D spectra of the Moyal decomposed 

shear stress as defined by Eq. (4.6.5). 
V1V2(K) indicates V1W2(K) , WXV2(K) 

W1W2(K) are the contributions to 
WAVE NO . indicates the magnitude of the 
wavenumber vector. Letters in the figure 
numbers correspond to the individual 
shear-flow simulations. 
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Fig. 4.14a. Contributing terms in the Reynolds 

stress equation defined as in Eq. (4.7.1) 
vs. St. Letters in the figure numbers cor- 
respond to the individual shear-flow simu- 
lations . 
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RUU STRESS EQUATION 


RUU STRESS EQUATION 



S*TIME S*TIME 


Fig. 4.14b. Contributing terms in the Reynolds Fig. 4.14c. Contributing terms in the R^^^ Reynolds 

stress equation defined as in Eq. (4.7.1) stress equation defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the individual shear -flow simulations. 
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Fig. 4.14d. Contributing terms in the R-^ Reynolds Fig. 4.14f. Contributing terms in the R^^ Reynolds 

stress equations defined as in Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters In the figure numbers correspond to the individual shear-flow simulations. 
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Flgi 4il4g. Contributing terms in tbe Reynolds Fig. 4.14h. Contributing terms in the R^^^ Re 3 molds 

stress equations defined as in Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the individual shear -flow simulations. 
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Fig. 4.15a. Contributing terms in the R-o Reynolds Fig. 4.15b. Contributing terms in the R ^2 Reynolds 

stress equations defined as xn Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the individual shear-flow simulations. 
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RVV STRESS EQUATION 


RVV STRESS EQUATION 



Fig. 4.15c. Contributing terms in the R,, Reynolds Fig. 4.15d. Contributing terms in the R Reynolds 

stress equations defined as xn Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the Individual shear-flow simulations. 




148 


RVV STRESS EQUATION 


RVV STRESS EQUATION 
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Fig. 4.15f. Contributing terms in the R „ Reynolds Fig. 4.15g. Contributing terms in the R ^2 Reynolds 

stress equations defined as xn Eq. (4,7. 1) stress equations defined as in Eq. (4. 7*1) 

vs. St. St. 

Letters in the figure nximbers correspond to the individual shear-flow simulations . 
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RVV STRESS EQUATION 


RWW STRESS EQUATION 



Fig. 4.15h. Contributing terms in the R„^ Reynolds Fig. 4.16a. Contributing terms in the R ^ Reynolds 
stress equations defined as in Eq. (4.7.1) stress equations defined as In Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the individual shear-flow simulations. 
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Fig. 4.16b. Contributing terms in the Reynolds Fig. 4.16c. Contributing terms in the R-^ Reynolds 

stress equations defined as in Eq . (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the individual shear— flow simulations. 
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Fig. 4.16d. Contributing terms in the Reynolds Fig. 4.16f. Contributing terms in the R^^ Reynolds 

stress equations defined as in Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the individual shear-flow simulations. 
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Fig. 4.16g. Contributing terms in the Reynolds Fig. 4.16h. Contributing terms in the R,^ Re 3 molds 

stress equations defined as in Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the individual shear -flow simulations. 
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RUV STRESS EQUATION 


RUV STRESS EQUATION 



Fig. 4.17a. Contributing terms in the Re 3 molds Fig. 4.17b. Contributing terms in the Reynolds 

stress equations defined as in Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 


Letters in the figure numbers correspond to the individual shear-flow simulations. 
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RUV STRESS EQUATION 



Fig. 4.17c. Contributing terms in the R ^2 Reynolds Fig. 4.17d. 
stress equations defined as in Eq. (4.7.1) 
vs. St. 


Contributing terms in the R ^2 Reynolds 
stress equations defined as in Eq. (4.7.1) 
vs. St. 


Letters in the figure numbers correspond to the individual shear-flow simulations. 
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RUV STRESS EQUATION 


RUV STRESS EQUATION 



S*TIME S*TIME 


Fig. 4.17f. Contributing terms in the 1^2 Reynolds Fig. 4.17g. Contributing terms in the ^2 Reynolds 
stress equations defined as in Eq. (4.7.1) stress equations defined as in Eq. (4.7.1) 

vs. St. vs. St. 

Letters in the figure numbers correspond to the Individual shear-flow simulations. 
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Fig. 4.17h. Contributing terms in the R ^2 Reynolds Fig. 4.18a. Contributing terms in the turbulent kinetic 
stress equations defined as in Eq. (4.7.1) energy equation vs. St. 

vs. St. 

Letters in the figure numbers correspond to the individual shear -flow simulations. 
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Fig. 4.18b. Contributing terms in the turbulent Fig. 4.18c. Contributing terms in the turbulent 

kinetic energy equation vs. St. kinetic energy equation vs. St. 

Letters in the figure numbers correspond to the individual shear-flow simulations. 
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Fig. 4.18d. Contributing terms in the turbulent Fig. 4.18f. Contributing terms in the turbulent 

kinetic energy equation vs. St. kinetic energy equation vs. St. 

Letters in the figure numbers correspond to the individual shear-flow simulations. 
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Fig. 4.19. Dissipation anisotropies vs. Reynolds 
stress anisotropies, all four nonzero 
components. 


Fig. 4.20. Second invariant of the compressible 
pressure-strain tensor divided by the 
sum of the second invariants of the Rotta 
and Fast pressure-strain tensors vs. the 
estimate of the least-squares fitted 
function. 




ROTTA MODEL (CONSTANT(C(1)) ROTTA TERM ALONE 


ROTTA MODEL (CONSTANT(C(1)) 2*EPS*DIJ INCLUDED 


CONSTANT/LEAST SQUARES EIT 
VS. VALUE OF THE CONSTANT 
o SHIFTED ORIGION ON THE ORDINATE 
8 - 


LEGEND 

o = terms in r(i,i) equation 
0 = terms in R(2.2) EQUATION 
A = terms in R(3.3) equation 
♦ = TERMS IN R(1.2) EQUATION 


CONSTANT/LEAST SQUARES FIT 
VS. VALUE OF THE CONSTANT 
SHIFTED ORIGION ON THE ORDINATE 


= TERMS IN Rf1,1) EQUATION 
= TERMS IN R(2.2) EQUATION 
"TERMS IN R(3,3) EQUATION 
"TERMS IN R(1,2) EQUATION 


® °o«o °o o 


A A a 

a'' 


2.5 3.0 3.5 

constant C(1) 


.5 4.0 

CONSTANT C(1) 


Fig. 4.21. Constant C.. in the Rotta pressure- 
strain model as used to replace just 
the Rotta pressure-strain term. 


Fig. 4.22. Constant C. in the Rotta pressure- 
strain model as used to replace just 
the Rotta pressure-strain term. 


An estimate of C^ from the least-squares fitted function is divided by C^ and plotted on the ordi- 
nate (the origin of the results for each of the four nonzero equations is shifted by two for clarity) . 


The value of 


is plotted on the abscissa. 
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Fig. 4.23. Constant A, in the fourth-rank tensor Fig. 4.24. Constant A„ in the Gibson-Launder 

model for the Fast pressure-strain terms. model for the Fast pressure-strain terms. 

An estimate of A from the least-squares fitted function is divided by A and plotted on the ordinate 
(the origin of the results for each of the four nonzero equations is shifted by two for clarity) . The 
value of A is plotted on the abscissa. 




RATIO OF PH i(3)11/PH 1(2)11 PSTRAIN VS. SL(11)/C 
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Fig. 4.25a. Ratio of compressible pressure-strain 
terms to Fast pressure-strain terms, 
3'^11^2^11 Mach number SL/c. 


Fig. 4.25b. Ratio of compressible pressure-strain 
terms to Fast pressure-strain terms, 
3 ^ 22 ^ 2^22 Mach number SL/c, 
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Fig. 4.25c. Ratio of compressible pressure-strain 
terms to Fast pressure-strain terms, 
3'^33^2*^33 Mach number SL/c. 


Fig. 4.25d. Ratio of compressible pressure-strain 
terms to Fast pressure-strain terms, 
3 ^ 12 ^ 2 '^H Mach number SL/c. 





FOURTH RANK TENSOR MODEL ON COMPRESSIBLE TERMS 


FOURTH RANK TENSOR MODEL ON FAST + COMP. TERMS 


constant/least squares fit 
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Fig. 4.26, Constant A., in the fourth-rank tensor Fig. 4.27. Constant k, in the fourth-rank tensor 

model used for the compressible pressure- model for trie sum of the compressible 

strain terms. and fast pressure-strain terms. 

An estimate of A^ from the least-squares fitted function is divided by A^ and plotted on the ordinate 
(the origin of the results for each of the four nonzero equations is shifted by two for clarity) . The 
value of A, is plotted on the abscissa. 


4.27. 
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Fig. 4.28. Test of structural similarity. The 

pressure-strain term is divided by the 
sum of all other contributing terms in 
the dynamic equations for the 

Reynolds stress anisotropy, ihis ratio 
is fit by least squares. The least- 
squares estimate is divided by the actual 
ratio and plotted on the ordinate. The 
ratio is plotted on the abscissa. 


Fig. 4.29. Constant in the structural similar- 

ity pressure-strain model, Eq . (4.12.5) 
for the total pressure-strain tensor. 

An estimate of from the least-squares 

fitted function is divided by and 

plotted on the ordinate (the origin of the 
results for each of the four nonzero equa- 
tions is shifted by two for clalrty) . The 
value of is plotted on the abscissa. 





APPENDIX 


TABULATION OF REDUCED DATA FROM 
THE SEVEN SHEAR FLOW SIMULATIONS 


All quantities are nondimensionalized on four constants; 


Density 


Sound Speed 


Reference length L = computational box side/2iT 

o 


Molecular viscosity p (constant) 


For example, the pressure is nondimensionalized on 


2 

P c 
o o 


Quantitative values for these dimensional quantities appear nowhere 
in the simulation because we always work with the nondimensional 
quantities. They have simply been used as guides to give the 
simulations Reynolds, Mach, and Shear numbers that are similar to 
those found in experiments. When these numbers are formed from the 
following length scales, etc., the dimensional quantities p , c , 
L^ drop out entirely, ° 

In this appendix, we use < > to indicate averages over the entire 
mesh . 
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M = < 


(u.Ui) 


1/2 


S = 


L 

dU o 
dy c 


Re. 


-1 


0 c T 
0^0 
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A 

4 
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00 
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A 

5 
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B 
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1.50E 

00 
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B 

6 

1.47E-01 

1.50E 

00 

2.00E-03 

C 

4 

3.12E-01 

3.00E 

00 
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C 

5 
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00 
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C 

6 
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D 

4 
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3.00E 
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D 

5 
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D 
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00 
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F 
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00 
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00 
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G 
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00 
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G 
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00 
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00 
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00 
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3.00E 

00 
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Taylor Integral Turbulent 

Micro Length Velocity 

Scale Scale 


F 1 ow 
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^ >,(!/« 




A 

4 
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00 
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00 
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00 
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1.25E 

00 
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3.64E-01 
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D 

4 
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00 

2.79E-01 

D 

5 

3.45E-01 
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00 
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0 
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1.50E 
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F 
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00 
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F 

5 
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F 
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00 
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G 
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3. 11E-01 

1.25E 

00 
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G 
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00 
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6 
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00 
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H 
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3. 18E-01 

1.43E 

00 
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Flow St 


SL/q Re <p’^/p2> 

o 
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01 
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4 
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00 

3.01E 

01 

7.21E-04 
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01 
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00 
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01 
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C 
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00 
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01 
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C 

6 
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01 
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01 
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D 

4 
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01 
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01 
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0 

5 
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01 
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01 
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D 

6 
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01 
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01 
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F 

4 
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01 
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01 
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F 

5 
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01 

6.63E 

01 
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F 

6 
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01 
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01 
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G 

4 

9.97E 

00 

2.93E 

01 

8.89E-03 

G 

5 

1.14E 

01 

3.34E 

01 

8.47E-03 

G 

6 

1.33E 

01 

3.79E 

01 

8.85E-03 

H 

4 

1.07E 

01 
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01 

4.81E-03 

H 

5 

1 .09E 

01 
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02 

5.53E-03 

H 

6 

1. 14E 

01 

1.21E 

02 

7.17E-03 


170 



asa:3:ooo'nTnoooooowa>0B>>> 


Decomposed Pressure Strain Tensor. 
Rotta Term. 
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Decomposed Pressure Strain Tensor. 
Past Term. 
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H 

4 

-1, 

.49E-02 

-4. 

23E-03 
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H 
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1 
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H 

6 
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1 
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6. 

94E-03 
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Decomposed Pressure Strain Tensor 


Compressible Term. 
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3:3cacooo*n'n-nooc?r>ocja>03W>>> 


Decomposed Reynolds Stress Tensor. 
Incompressible Part. 
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4.81E-03 
4.27E-03 
2.S7E-02 
2.27E-02 
2.04E-02 
8.54E-03 
7.86E-03 
7.21E-03 
7.74E-03 
6.87E-03 
6.22E-03 
1.87E-02 
1.54E-02 
1 .30E-02 
1.34E-02 
1.38E-02 
1.43E-02 


1.91E-03 
1.83E-03 
1.81E-03 
1 .09E-02 
1.07E-02 
1.08E-02 
4.92E-02 
4.84E-02 
4.81E-02 
1.90E-02 
1.90E-02 
1 .94E-02 
1.96E-02 
1.88E-02 
1 .85E-02 
3.77E-02 
3.57E-02 
3.42E-02 
2.74E-02 
2.93E-02 
3. 19E-02 
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3:3:3:ooo-n-n-naocjooocoa300>>> 


Decomposed Reynolds Stress Tensor 
Compressible Part. 


R 


D 


2 

p^c 
0 o 


Flow 


St 


i=1>j=2 i=1.j=1 i=2,j=2 i=3,j=3 




4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 


-1.30E-06 
-2.27E-05 
-3. 18E-05 
-2.40E-05 
-1, 14E-04 
-2. 13E-04 
-1.57E-03 
-2.69E-03 
-2.02E-03 
-4.88E-04 
-9.65E-04 
-7. 13E-04 
-1.36E-03 
-1. 16E-03 
-1.25E-03 
-1.39E-03 
-2.30E-03 
-1.67E-03 
-5.57E-04 
-1. 17E-03 
-8.89E-04 


4.65E-05 
4.69E-05 
7. 17E-05 
3.21E-04 
2.66E-04 
3.93E-04 
3.97E-03 
4.69E-03 
3.29E-03 
1.21E-03 
1.61E-03 
1.05E-03 
2.56E-03 
2.54E-03 
2.43E-03 
3.38E-03 
3.93E-03 
2.52E-03 
1.44E-03 
2.01E-03 
1.35E-03 


-3.08E-05 
7.29E-06 
-6.68E-07 
-1.90E-04 
4.54E-05 
4.23E-05 
1.13E-03 
2.48E-03 
4.65E-03 
2.86E-04 
9.38E-04 
I.73E-03 
2.27E-03 
2.55E-03 
3.32E-03 
9. 19E-04 
2.08E-03 
3.89E-03 
3.97E-04 
1. 17E-03 
2.23E-03 


2.73E-05 
6.52E-06 
3.79E-05 
1.52E-04 
5.65E-05 
2.88E-04 
1.38E-03 
1 .80E-03 
1.90E-03 
4.73E-04 
4.84E-04 
5.42E-04 
1.08E-03 
1.44E-03 
9.04E-04 
1.22E-03 
1.50E-03 
1.50E-03 
5.07E-04 
5.53E-04 
6.68E-04 


175 



Homogeneous Dissipation Tensor. 



Re, ^ < u . , u , , > 
b i,k j,k 



2 


c 

o 


Flou St i=1,j=1 i=2,j=2 i=3,j=3 

***»********»*«****«»*»*****«»»*»**»*»«»**»*»*»»*»**»»*»***» 


A 

4 

2.09E-03 

2.59E-04 

1. 16E-03 

A 

5 

2. 14E-03 

2.05E-04 

9.69E-04 

A 

6 

2. 14E-03 

1.74E-04 

8.39E-04 

B 

4 

1 .01E-02 

2.45E-03 

6.41E-03 

B 

5 

9.42E-03 

1 .9SE-03 

5.53E-03 

B 

6 

9.08E-03 

1.74E-03 

4.92E-03 

C 

4 

8.47E-02 

2.53E-02 

5.26E-02 

C 

5 

8.09E-02 

2. 19E-02 

4.65E-02 

C 

6 

7.86E-02 

2.09E-02 

4.20E-02 

0 

4 

3. 13E-02 

5.68E-03 

1.45E-02 

0 

5 

3.34E-02 

5.38E-03 

1.32E-02 

D 

6 

3.50E-02 

5.46E-03 

1.27E-02 

F 

4 

5.26E-02 

7.93E-03 

1 .84E-02 

F 

5 

6.22E-02 

7. 13E-03 

1.66E-02 

F 

6 

7.02E-02 

7.32E-03 

1.54E-02 

G 

4 

7.90E-02 

2.06E-02 

4.77E-02 

G 

5 

7.32E-02 

1 .69E-02 

4.04E-02 

G 

6 

6.90E-02 

1.57E-02 

3.49E-02 

H 

4 

2.59E-02 

5.72E-03 

1 .20E-02 

H 

5 

2.98E-02 

6. 10E-03 

1.21E-02 

H 

6 

3.33E-02 

6.79E-03 

1.27E-02 
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3:a:a:c>CTo-Tn-noocJoooeoa>ro>>> 


Homogeneous Dissipation Tensor 






X 



2 

c 

o 


Flow 


St 


i=1»j=2 i=1,j=3 i=2,j=3 




4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 

4 

5 

6 


-4 

,54E- 

-04 

-1. 

,67E- 

■06 

8. 

39E- 

■06 

-3 

,81E- 

-04 

-3. 

21E- 

-06 

2. 

70E- 

-06 

-3 

,22E- 

-04 

-5. 

. 15E- 

■06 

-1. 

66E- 

■06 

-2 

,63E- 

-03 

3. 

, 13E- 

-05 

3. 

20E- 

■05 

-2 

, 12E- 

-03 

-8. 

,28E- 

■06 

3. 

20E- 

■05 

-1 

, 83E- 

-03 

-3. 

.21E- 

■05 

6. 

52E- 

■06 

-2 

• 28E- 

-02 

3. 

.07E- 

-04 

4. 

88E- 

-04 

-1 

90E- 

-02 

-1. 

26E- 

■04 

3. 

18E- 

-04 

-1 

. 64E- 

-02 

-5. 

95E- 

-04 

4. 

77E- 

■05 

-7 

,82E- 

-03 

4. 

,96E- 

■06 

9. 

38E- 

-05 

-7 

• 25E- 

-03 

1, 

,75E- 

-04 

-7. 

21E- 

-06 

-6 

,56E- 

-03 

5. 

49E- 

-06 

-1 . 

82E- 

-04 

-1 

,21E- 

-02 

2. 

,04E- 

-05 

1. 

40E- 

-04 

-1 

, 19E- 

-02 

3. 

68E- 

-04 

-2. 

50E- 

-05 

-1 

. 14E- 

-02 

2, 

,07E- 

■04 

-1. 

83E- 

■04 

-2 

. 12E- 

■02 

2. 

,03E- 

-04 

5. 

OOE- 

-04 

-1 

.71E- 

-02 

-1. 

77E- 

-04 

2. 

76E- 

■04 

-1 

• 42E- 

-02 

-5. 

, 15E- 

-04 

-1. 

11E- 

-04 

-6 

,68E- 

-03 

9. 

.38E- 

-06 

8. 

32E- 

■05 

-6 

• 79E- 

-03 

1, 

87E- 

-04 

-1. 

50E- 

■05 

-6 

• 64E- 

-03 

8. 

■ 01E- 

-05 

-1. 

58E- 

-04 
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Intensities of the decomposed pressure fields. 


Rotta Fast Compr. 


F 1 ow 






A 

4 

5.95E-06 

2.04E-05 

7.71E-05 

A 

5 

5.53E-06 

2.45E-05 

1 .21E-04 

A 

6 

5.04E-06 

2.24E-05 

1.68E-04 

B 

4 

1.75E-04 

1 . 58E-04 

6. 33E-04 

B 

5 

1.62E-04 

1 .82E-04 

1.02E-03 

B 

6 

1.57E-04 

1.53E-04 

1 .22E-03 

C 

4 

4.08E-03 

1 .80E-03 

1.59E-02 

c 

5 

3.97E-03 

2.00E-03 

1.67E-02 

c 

6 

3.93E-03 

1.7SE-03 

2. 12E-02 

D 

4 

7. 10E-04 

5. 11E-04 

4. 16E-03 

D 

5 

7.36E-04 

6. 18E-04 

4.65E-03 

D 

6 

7.59E-04 

5.95E-04 

6.83E-03 

F 

4 

1.20E-03 

1. 17E-03 

1 .03E-02 

F 

5 

1 .24E-03 

1.32E-03 

1.05E-02 

F 

6 

1 .23E-03 

1. 16E-03 

1. 13E-02 

G 

4 

2. 55E-03 

1 . 51E-03 

1 -29E-02 

G 

5 

2.32E-03 

1.63E-03 

1.32E-02 

6 

6 

2.20E-03 

1.41E-03 

1.69E-02 

H 

4 

1.41E-03 

6.33E-04 

5.46E-03 

H 

5 

1.64E-03 

8.05E-04 

6.29E-03 

H 

6 

1.90E-03 

8. 13E-04 

9.46E-03 
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